Welcome, guest | Sign In | My Account | Store | Cart

If you are calling factorial() many times, it is worth to cache the computed factorials.

Python, 29 lines
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
_FAC_TABLE = [1, 1]
def factorial(n):
    if n < len(_FAC_TABLE):
        return _FAC_TABLE[n]

    last = len(_FAC_TABLE) - 1
    total = _FAC_TABLE[last]
    for i in range(last + 1, n + 1):
        total *= i
        _FAC_TABLE.append(total)

    return total


def main():
    from timeit import Timer
    print("factorial:      %.5f s" %
            Timer("factorial(500)", "from __main__ import factorial")
            .timeit(1000))

    import math
    if hasattr(math, "factorial"):
        print("math.factorial: %.5f s" %
                Timer("factorial(500)", "from math import factorial")
                .timeit(1000))


if __name__ == "__main__":
    main()

My timeit results:

factorial:      0.00206 s
math.factorial: 0.35708 s

When computing factorial(n), any factorial(i) for i < n will be also precomputed. And factorial(n) will be reused when computing factorial(n + 1) and higher.

6 comments

Marc Liddell 13 years, 10 months ago  # | flag

All you're testing is how fast python can do the first calculation and then return an element in a list 999 times.

By that same logic instead of the complicated process you're using you may as well combine both ideas to make the far simpler:

def mfactorial(n, _factcache={}):
  if n not in _factcache:
    _factcache[n] = math.factorial(n)
  return _factcache[n]

factorial:      0.00080 s
mfactorial:      0.00059 s
math.factorial: 0.25890 s
Ivo Danihelka (author) 13 years, 10 months ago  # | flag

You are right that the test does not measure all usages. Your mfactorial() has a disadvantage. It does not store the computed factorial(1) .. factorial(n-1).

Daniel Stutzbach 13 years, 10 months ago  # | flag

For what it's worth, Python 3.2 will include a much faster factorial function than earlier versions of Python (see Issue8692).

David Lambert 13 years, 9 months ago  # | flag

To compute factorials faster for large arguments, use Stirling's Approximation (http://en.wikipedia.org/wiki/Stirling%27s_approximation):

fN = float(N)
Log(N!) = N*math.log(fN) + math.log(2*math.pi*fN)/2 - fN + \
      (fN**-1)/12 - (fN**-3)/360 + (fN**-5)/1260 - (fN**-7)/1680 + (fN**-9)/1188

As N grows, the negative powers of N improve the approximation.
This approximation is so good that:

  • For N > 11, the (fN^-9) term hits float's limit of precision.
  • For N > 17, the (fN^-7) term hits float's limit of precision.
  • For N > 32, the (fN^-5) term hits float's limit of precision.
  • For N around 94, the (fN^-3) term hits float's limit of precision.
  • For N around 670, the (fN^-1) term hits float's limit of precision.

When a term hits float's limit of precision, you can omit it from the approximation.

This program shows the accuracy of Stirling's Approximation:

import math

# Create output file.
sFile = 'C:\Documents\Desktop\Factorial_Float.tsv'
try:
    oFile = open(sFile, 'w')
except:
    print "Problem opening file " + sFile + ".\nGoodbye"
    quit()

# Column headers.
sOutput = "N\tExact\t0\t1\t2\t3\t4\t5\n"
oFile.write( sOutput )

# Calculate exact & approx values of Ln( N! ) for N = 1 to 1000.
# Long variables have prefix "i", float variables have prefix "f".

# Exact value of Ln( N! ).
fExact = 0
for iN in range(1,1000):
    # Float version of iN.
    fN = float(iN)
    fExact = fExact + math.log(fN)
    sOutput = str(iN) + "\t" + str(fExact)

    # f1 = 1st term of Stirling's Approx.
    f1 = math.log( 2*math.pi*fN )/2 + fN * math.log( fN ) - fN
    # fDiff always = difference between exact & approx values of Ln( N! )
    fDiff = f1-fExact
    sOutput = sOutput + "\t" + str(abs(fDiff))

    # f2 = 2nd term of Stirling's Approx.
    # dDelta always = next term in Stirling's Approx.
    fDelta = (fN**-1)/12
    f2 = f1 + fDelta
    fDiff = fDiff + fDelta
    sOutput = sOutput + "\t" + str(abs(fDiff))

    # f3 = 3rd term of Stirling's Approx.
    fDelta = -(fN**-3)/360
    f3 = f2 + fDelta
    fDiff = fDiff + fDelta
    sOutput = sOutput + "\t" + str(abs(fDiff))

    # f4 = 4th term of Stirling's Approx.
    fDelta = (fN**-5)/1260
    f4 = f3 + fDelta
    fDiff = fDiff + fDelta
    sOutput = sOutput + "\t" + str(abs(fDiff))

    # f5 = 5th term of Stirling's Approx.
    fDelta = -(fN**-7)/1680
    f5 = f4 + fDelta
    fDiff = fDiff + fDelta
    sOutput = sOutput + "\t" + str(abs(fDiff))

    # f6 = 6th term of Stirling's Approx.
    fDelta = (fN**-9)/1188
    f6 = f5 + fDelta
    fDiff = fDiff + fDelta
    sOutput = sOutput + "\t" + str(abs(fDiff))

    oFile.write( sOutput + "\n" )

oFile.close()

print "Goodbye"
Giannis Fysakis 13 years, 9 months ago  # | flag

David Lambert this method seems pretty accurate ,...but when i give value to N greater than 170 i get an overflowError

Example code:

import math

for N in xrange(1,10):
    fN=float(N)
    FactorialOfN=N*math.log(fN) + math.log(2*math.pi*fN)/2 - fN + (fN**-1)/12 \
    - (fN**-3)/360 + (fN**-5)/1260 - (fN**-7)/1680 + (fN**-9)/1188
    Result = math.e **FactorialOfN
    print "Factorial of %d" %N ," = " , Result


for N in xrange(10,190,10):
    fN=float(N)
    FactorialOfN=N*math.log(fN) + math.log(2*math.pi*fN)/2 - fN + (fN**-1)/12 \
    - (fN**-3)/360 + (fN**-5)/1260 - (fN**-7)/1680 + (fN**-9)/1188
    Result = math.e **FactorialOfN
    print "Factorial of %d" %N ," = " , Result

Example Error:

~$ ./Stirlingfactorial.py 
    Factorial of 1  =  1.00053439504
    Factorial of 2  =  2.00000108765
    Factorial of 3  =  6.00000004849
    Factorial of 4  =  24.0000000092
    Factorial of 5  =  120.000000004
    Factorial of 6  =  720.000000003
    Factorial of 7  =  5040.0
    Factorial of 8  =  40320.0
    Factorial of 9  =  362880.0
    Factorial of 10  =  3628800.0
    Factorial of 20  =  2.43290200818e+18
    Factorial of 30  =  2.65252859812e+32
    Factorial of 40  =  8.15915283248e+47
    Factorial of 50  =  3.04140932017e+64
    Factorial of 60  =  8.32098711274e+81
    Factorial of 70  =  1.197857167e+100
    Factorial of 80  =  7.15694570463e+118
    Factorial of 90  =  1.48571596448e+138
    Factorial of 100  =  9.33262154439e+157
    Factorial of 110  =  1.58824554152e+178
    Factorial of 120  =  6.68950291345e+198
    Factorial of 130  =  6.46685548922e+219
    Factorial of 140  =  1.34620124757e+241
    Factorial of 150  =  5.71338395644e+262
    Factorial of 160  =  4.71472363599e+284
    Factorial of 170  =  7.25741561531e+306
    Traceback (most recent call last):
      File "./Stirlingfactorial.py", line 16, in <module>
        Result = math.e **FactorialOfN
    OverflowError: (34, 'Numerical result out of range')

Can you suggest a way to get rid of this error :) ?

Thanks in advance :-)

Giannis Fysakis 13 years, 9 months ago  # | flag

i understand it's a really big number for the Computer... to represent it ... so i assume we can keep it in it's "log" representation.... for "mathematistic" manipulations and "stuff" ...

Created by Ivo Danihelka on Sun, 23 May 2010 (MIT)
Python recipes (4591)
Ivo Danihelka's recipes (2)

Required Modules

  • (none specified)

Other Information and Tasks