Problem with extremely small real number

[Arnaud Delobelle]

On Sep 3, 3:56 pm, Andrea <aciru... at gmail.com> wrote:
> Hi,
> I need to calculate this probability P!/{n \choose p}, varying both n
> and p, n takes values in this range [512:1024] and p in [2:12].
> So i write this code in python:
>
> def factorial(n):
>         result=1
>         if n==0: return 1
>         for i in xrange(1, abs(n)+1):
>                 result = i*result
>         if n >= 0:
>                 return result
>         else:
>                 return -result
> def binomial(n, k):
>         assert n>0 and isinstance(n, (int, long)) and isinstance(k,
> (int,long))
>         if k < 0 or k > n:
>                 return 0
>         if k == 0 or k == n:
>                 return 1
>         return factorial(n) // (factorial(k) * factorial(n-k))
>
> I want to call factorial(2)//binomial(1024,2) for example, in this way
> trivially I obtain 0 as probability, how can I obtain the probability
> forcing this division to output extremely small real numbers????
>
> I want to produce a program that provide a set of data of this
> probability on varying n and p, in order to plot a surface with this
> data set with gnuplot, any comments or suggestion?
> thanks in advance,
> Andrea

You're using integer division, that's why you probabilities are 0!
In python (up to 2.5 at list), if one integer/long is divided by
another, integer division is used.
But if one of the operands is a float, then the other one is converted
to a float and float division is used.

* convert your integers to floats first.
* Don't use //, use /
(// is 'true division', it will always give you the rounded answer)

[This can cause confusion, and I think it has been changed from v2.6
on - which I haven't got to try so correct me if  I'm wrong]

Thus:
\
>>> float(factorial(2)) / binomial(1024, 2)
3.8184261974584555e-06

BTW
p! / nCp = p! / (n! / p!(n-p)!) = (p!)^2 * (n! / (n-p)!) = (p!)^2 /
\prod_{i=n-p+1}^n i

Therefore:

def product(a, b):
    return reduce(lambda x, y: x*y, xrange(a, b+1))

def f(n, p):
    return float(product(1, p))**2 / product(n-p+1, n)

>>> f(1024, 2)
3.8184261974584555e-06

This will be more efficient.

--
Arnaud