# Calculating factorial

Tom Loredo loredo at spacenet.tn.cornell.edu
Thu Dec 7 23:00:23 CET 2000

```Howdy-

The functions below are Python versions of functions appearing in the
*Numerical Recipes* books by Press, Teukolsky, etc..  factrl(n)
returns n! as a float; factln(n) returns the natural logarithm
of n! as a float.  Both store values calculated for lowish n to
speed subsequent calls.  They are useful for calculations in
statistics, but probably not what you'd want to use for combinatorics.
Just wrote them a couple days ago; use at your own risk!

This is pure Python for portability; for speed you might consider
using gamma function calls from the cephes library.

Peace,
Tom Loredo

import Numeric as N
from math import log, exp

gammln_cof = N.array([76.18009173, -86.50532033, 24.01409822,
-1.231739516e0, 0.120858003e-2, -0.536382e-5])
gammln_stp = 2.50662827465

def gammln(xx):
"""Logarithm of the gamma function."""
global gammln_cof, gammln_stp
x = xx - 1.
tmp = x + 5.5
tmp = (x + 0.5)*log(tmp) - tmp
ser = 1.
for j in range(6):
x = x + 1.
ser = ser + gammln_cof[j]/x
return tmp + log(gammln_stp*ser)

def factrl(n, ntop=0, prev=N.ones((33),N.Float)):
"""Factorial of n.
The first 33 values are stored as they are calculated to
speed up subsequent calculations."""
if n < 0:
raise ValueError, 'Negative argument!'
elif n <= ntop:
return prev[n]
elif n <= 32:
for j in range(ntop+1, n+1):
prev[j] = j * prev[j-1]
ntop = n
return prev[n]
else:
return exp(gammln(n+1.))

def factln(n, prev=N.array(101*(-1.,))):
"""Log factorial of n.
Values for n=0 to 100 are stored as they are calculated to
speed up subsequent calls."""
if n < 0:
raise ValueError, 'Negative argument!'
elif n <= 100:
if prev[n] < 0:
prev[n] = gammln(n+1.)
return prev[n]
else:
return gammln(n+1.)

```