44

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.

Martijn Pieters
  • 1,048,767
  • 296
  • 4,058
  • 3,343
user2312186
  • 453
  • 1
  • 4
  • 8

5 Answers5

43

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)
Martijn Pieters
  • 1,048,767
  • 296
  • 4,058
  • 3,343
  • What is the large number class in Python called? because `170!` is well above the limits of a 64-bit integer. I assume the type is coerced to that? If so, why is the `L` still appended to the end? – Hunter McMillen Apr 23 '13 at 16:31
  • 3
    @HunterMcMillen: Python long integers are only limited by OS memory allocation. The `L` signifies we moved beyond the maximum C integer size for the platform (64-bit in my case) to show long integers are being used; the transition is automatic. See http://docs.python.org/2/library/stdtypes.html#numeric-types-int-float-long-complex – Martijn Pieters Apr 23 '13 at 16:33
  • Ahh, I see. Most other languages have `BigInteger` or something similar. Thanks. – Hunter McMillen Apr 23 '13 at 16:36
7

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.

thepith
  • 67
  • 6
Pyrce
  • 8,296
  • 3
  • 31
  • 46
4

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.

xylon97
  • 264
  • 1
  • 5
  • Hi guys I tried to use decimal or scipy but i got ImportError: No module named scipy.stats the same error for decimal – user2312186 Apr 23 '13 at 17:08
  • what version are you running? `decimal` library is supposed to be standard for all. perhaps you have deleted it. re-install python and it will fix it. as for scipy, it is a third-party package, so you'll have to download it from scipy.org – xylon97 Apr 23 '13 at 17:12
  • I am going to use poisson like this: depart = depart + Poisson(3) so that inter-arrival times are selected by Poisson. but it does not work. – user2312186 Apr 23 '13 at 17:45
2

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 ]
Charles Brunet
  • 21,797
  • 24
  • 83
  • 124
  • +1 for `scipy.stats`, but the factorial function scipy offers will still `inf` out at large argument. – DSM Apr 23 '13 at 16:46
  • Thanks. I am going to generate Poisson numbers as generating packet time in the network. How can i use this function for generating a Poisson number each continually? – user2312186 Apr 23 '13 at 16:52
0

If you get this error trying to do division with very large numbers, use the floor division operator(//) instead of /.

/ always returns a float whether or not the quotient is a whole number, but // will always return an int because it ignores any remainder.

Floor division works well when you only need the nearest whole number, or when you know that there will be no remainder.

For example:

n = 5479174226627386185132548618631889901184058567721039026663537776859958572142823922154018004124214331284656070023110969954188349772902119804704702598117413535893235324950786546151906304553125886840777572957470085282511148870434624033780058878942475280220905204940912590689793763351716319910433886468914385951305420981218341758360474092776301397873143398400716102216234362319599979129849103060426402698870872632043715074191472826989513693244855131547410762874345810755472991691149113189232066769661520546071424890254793363936991126744220856839621014660137423821784120362098968560498064575399109509809144075148032776900847515603985610957107880554355044275212955190073519527570422912054125513731742317779855206190407418289517246109264745652852110081421227227070961742013273408514723195685741381697286381613758857699193539587512284537797999746347217299424379359513750952221691770954364716137930556670199833832758082670990435174217159016818867385870452581971005123842835195274660590963683055276508447399891075634175690730834545095637653926683291205950640671169613344094321183481673015248807254721861961606825265236717120630995029306863300152797740979622723077992597531418730154060586722026762649940060327662106126471869300626956821351541276596745066988430327015357076633723108282682951496900104227774031265608955533139148284989331915052933624499658651051440828202238489960326767828269045239036246524560280599359429362948869990238325225414863546810033093049656312737083470849175724716280988160656752494507880308537719866213986281428744395701992495890900685006572556038647427947211000473774809669487947261242971504873560862774351759068398609198938063648190986411623475678601785762035198084851677024720418537109078696439281976199801473556400242587983954277326038907443266131383984144369164421859162116568291430549701757135817698503903004323670868472322098027272184621482168087031650528289698350186047041813000034365448608048428773660686833199312751498967864183755218605831200984734033178462911134687785324823763323028015552729641117736966587793077203826296036041919194725026292750532442085716784596992382527663553269382473454979408093249648685286746030229491081455014614036595064858640568694698486979304728018771546672574024918565616067368338166156732338857331627951808169957419646902230036338741999953636036848632388844148127317892455668907319296843822473855127548875966047711681567056047142200887652753813947766711904166025820768269034575968052119510503904027778183601434729384390496315574344178298921662230866717650218116213036528198477594125535378602905735898504877520174927689474709976093921901862335223731814101741854177482601528306067859536805771591038922076905431603327032115233190666894167725361120047719134862842475859620392485715336420653304556079280724582370679687097844646877141209557539992623586632733333757970772758010549948610992979279029728363939092740579103618926669219471007734885307035386951609417262492554868268269616601555317177098129073353194374818770380635286890143608518924421556937379652833811709922933535682467614307592250619010586908268258080894148787717354387660022563306126433196292155158692150378675047252139601642814517370404428943577553004710999041554686557359791855943269834124640464693626652005988300338967392066239532981003255609078688925846590550467734155623629681262423331018717094741382157740779203008752310807956601087650996558858762789654321362720766176441413026718299435000087843048352098779842582683343460060362948117211949633698490613691410568825978953300376848779528540637469856488967807293836747161816004817653656673962008997810605012877379655141001116834025783685768201503435029018375437319930124960810208705813192640556535800129029066957647413650553915827395477280121732521449721923559779963923510646372389370353922598568377127428196541696965429553407721148989373822271488


divisor = 2

if n % divisor == 0:
    # We know the number is evenly divisible by the divisor, so we can use // without losing any possible remainder
    quotient = n // divisor  # This will work just fine

# This line will give an error since / always returns a float
quotient = n / divisor
# OverflowError: integer division result too large for a float
Captain Jack Sparrow
  • 971
  • 1
  • 11
  • 28