Problem with extremely small real number
Arnaud Delobelle
arnodel at googlemail.com
Mon Sep 3 17:52:01 CEST 2007
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
More information about the Python-list
mailing list