Operations for Long and Float in Python

Question

I'm trying to compute this:

from scipy import *
3600**3400 * (exp(-3600)) / factorial(3400)

the error: unsupported long and float

Answer 1

Well the error is coming about because you are trying to multiply

3600**3400

which is a long with

exp(-3600)

which is a float.

But regardless, the error you are receiving is disguising the true problem. It seems exp(-3600) is too big a number to fit in a float anyway. The python math library is fickle with large numbers, at best.

Answer 2

Try using logarithms instead of working with the numbers directly. Since none of your operations are addition or subtraction, you could do the whole thing in logarithm form and convert back at the end.

Answer 3

You could try using the Decimal object. Calculations will be slower but you won't have trouble with really small numbers.

from decimal import Decimal

I don't know how Decimal interacts with the scipy module, however.

This numpy discussion might be relevant.

Answer 4

Computing with numbers of such magnitude, you just can't use ordinary 64-bit-or-so floats, which is what Python's core runtime supports. Consider gmpy (do not get the sourceforge version, it's aeons out of date) -- with that, math, and some care...:

>>> e = gmpy.mpf(math.exp(1))
>>> gmpy.mpz(3600)**3400 * (e**(-3600)) / gmpy.fac(3400)
mpf('2.37929475533825366213e-5')

(I'm biased about gmpy, of course, since I originated and still participate in that project, but I'd never make strong claims about its floating point abilities... I've been using it mostly for integer stuff... still, it does make this computation possible!-).

Answer 5

exp(-3600) is too smale, factorial(3400) is too large:

In [1]: from scipy import exp

In [2]: exp(-3600)
Out[2]: 0.0
In [3]: from scipy import factorial

In [4]: factorial(3400)
Out[4]: array(1.#INF)

What about calculate it step by step as a workaround(and it makes sense to check the smallest and biggest intermediate result):

from math import exp
output = 1
smallest = 1e100
biggest = 0
for i,j in izip(xrange(1, 1701), xrange(3400, 1699, -1)):
    output = output * 3600 * exp(-3600/3400) / i
    output = output * 3600 * exp(-3600/3400) / j
    smallest = min(smallest, output)
    biggest = max(biggest, output)
print "output: ", output
print "smallest: ", smallest
print "biggest: ", biggest

output is:

output:  2.37929475534e-005
smallest:  2.37929475534e-005
biggest:  1.28724174494e+214