ISO module for binomial coefficients, etc.

[Daniel Stutzbach]

On Sat, Jan 23, 2010 at 4:55 PM, kj <no.email at please.post> wrote:

> Before I go off to re-invent a thoroughly invented wheel, I thought
> I'd ask around for some existing module for computing binomial
> coefficient, hypergeometric coefficients, and other factorial-based
> combinatorial indices.  I'm looking for something that can handle
> fairly large factorials (on the order of 10000!), using floating-point
> approximations as needed, and is smart about optimizations,
> memoizations, etc.
>

I don't have code for any of the other stuff, but I have some efficient code
(see below) for computing exact integer factorials.  It's a
divide-and-conquer algorithm.  I'm not sure of the exact time complexity,
but it appears to somewhere between n*log n and n**2.

python2.6 -m timeit -s 'from utils import factorial' 'x = factorial(10000)'
10 loops, best of 3: 27.7 msec per loop

If that's not fast enough, I suggest using the gamma function to compute
floating-point approximations, as gamma(n) == (n-1)!.  Gamma isn't yet
included in the Python standard library (http://bugs.python.org/issue3366),
but you can use the ctypes module to get it from the system C library on
many platforms.

For computing binomial coefficients, you can use the lgamma function (log of
gamma), also found in my system C libraries.  Since choose(n, k) =
exp(lgamma(n) - lgamma(k) - lgamma(n-k)).

If you go with floating point, watch out for floating point rounding error.
;-)

def factorial(n):
    "http://www.luschny.de/math/factorial/binarysplitfact.html"
    p = 1
    r = 1
    p, r = factorial_loop(n, p, r)
    return r << n_minus_num_of_bits(n)

def factorial_loop(n, p, r):
    if n <= 2: return p, r
    p, r = factorial_loop(n // 2, p, r)
    p *= part_product(n // 2 + 1 + ((n // 2) & 1),
                      n - 1 + (n & 1))
    r *= p
    return p, r

def part_product(n, m):
    if (m <= (n+1)): return n
    if (m == (n+2)): return n*m
    k = (n+m) // 2
    if ((k & 1) != 1):
        k -= 1
    return part_product(n, k) * part_product(k+2, m)

def n_minus_num_of_bits(v):
    w = v
    if w >= 2**32-1:
        raise OverflowError
    w -= (0xaaaaaaaa & w) >> 1
    w = (w & 0x33333333) + ((w >> 2) & 0x33333333)
    w = w + (w >> 4) & 0x0f0f0f0f
    w += w >> 8
    w += w >> 16
    return v - (w & 0xff)

--
Daniel Stutzbach, Ph.D.
President, Stutzbach Enterprises, LLC <http://stutzbachenterprises.com>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/python-list/attachments/20100127/50c355ce/attachment-0001.html>