I tried to calculate poisson distribution in python as below:

``````p = math.pow(3,idx)
depart = math.exp(-3) * p
depart = depart / math.factorial(idx)
``````

idx ranges from 0

But I got `OverflowError: long int too large to convert to float`

I tried to convert depart to `float` but no results.

Factorials get large real fast:

``````>>> math.factorial(170)
7257415615307998967396728211129263114716991681296451376543577798900561843401706157852350749242617459511490991237838520776666022565442753025328900773207510902400430280058295603966612599658257104398558294257568966313439612262571094946806711205568880457193340212661452800000000000000000000000000000000000000000L
``````

Note the `L`; the factorial of 170 is still convertable to a float:

``````>>> float(math.factorial(170))
7.257415615307999e+306
``````

but the next factorial is too large:

``````>>> float(math.factorial(171))
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
OverflowError: long int too large to convert to float
``````

You could use the `decimal` module; calculations will be slower, but the `Decimal()` class can handle factorials this size:

``````>>> from decimal import Decimal
>>> Decimal(math.factorial(171))
Decimal('1241018070217667823424840524103103992616605577501693185388951803611996075221691752992751978120487585576464959501670387052809889858690710767331242032218484364310473577889968548278290754541561964852153468318044293239598173696899657235903947616152278558180061176365108428800000000000000000000000000000000000000000')
``````

You’ll have to use `Decimal()` values throughout:

``````from decimal import *

with localcontext() as ctx:
ctx.prec = 32  # desired precision
p = ctx.power(3, idx)
depart = ctx.exp(-3) * p
depart /= math.factorial(idx)
``````

When `idx` gets large either the `math.pow` and/or the `math.factorial` will become insanely large and be unable to convert to a floating value (`idx=1000` triggers the error on my 64 bit machine). You’ll want to not use the math.pow function as it overflows earlier than the built in `**` operator because it tries to keep higher precision by float converting earlier. Additionally, you can wrap each function call in a `Decimal` object for higher precision.

Another approach when dealing with very large numbers is to work in the log scale. Take the log of every value (or calculate the log version of each value) and perform all required operations before taking the exponentiation of the results. This allows for your values to temporary leave the floating domain space while still accurately computing a final answer that lies within floating domain.

``````3 ** idx  =>  math.log(3) * idx
math.exp(-3) * p  =>  -3 + math.log(p)
math.factorial(idx)  =>  sum(math.log(ii) for ii in range(1, idx + 1))
...
math.exp(result)
``````

This stays in the log domain until the very end so your numbers can get very, very large before you’ll hit overflow problems.

Try using the decimal library. It claims to support arbitrary precision.
`from decimal import Decimal`

Also, you don’t need to use `math.pow`. `pow` is in-built.

The `scipy` module could help you.

scipy.misc.factorial is a factorial function that can use the gamma function approximation to calculate the factorial, and returns the result using floating points.

``````import numpy
from scipy.misc import factorial

i = numpy.arange(10)
print(numpy.exp(-3) * 3**i / factorial(i))
``````

Gives:

``````[ 0.04978707  0.14936121  0.22404181  0.22404181  0.16803136  0.10081881
0.05040941  0.02160403  0.00810151  0.0027005 ]
``````

There is also a module to calculate Poisson distributions. For example:

``````import numpy
from scipy.stats import poisson

i = numpy.arange(10)
p = poisson(3)
print(p.pmf(i))
``````

Gives:

``````[ 0.04978707  0.14936121  0.22404181  0.22404181  0.16803136  0.10081881
0.05040941  0.02160403  0.00810151  0.0027005 ]
``````