Calculating factorial

[Tom Loredo]

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.)