<div class="gmail_quote">On Sat, Jan 23, 2010 at 4:55 PM, kj <span dir="ltr"><no.email@please.post></span> wrote:<br><blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">

Before I go off to re-invent a thoroughly invented wheel, I thought<br>
I'd ask around for some existing module for computing binomial<br>
coefficient, hypergeometric coefficients, and other factorial-based<br>
combinatorial indices.  I'm looking for something that can handle<br>
fairly large factorials (on the order of 10000!), using floating-point<br>
approximations as needed, and is smart about optimizations,<br>
memoizations, etc.<br></blockquote><div><br>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.<br>
<br>python2.6 -m timeit -s 'from utils import factorial' 'x = factorial(10000)'<br>
10 loops, best of 3: 27.7 msec per loop<br>
<br>
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 (<a href="http://bugs.python.org/issue3366">http://bugs.python.org/issue3366</a>), but you can use the ctypes module to get it from the system C library on many platforms.<br>
<br>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)).<br><br>If you go with floating point, watch out for floating point rounding error. ;-)<br>
<br>def factorial(n):<br>    "<a href="http://www.luschny.de/math/factorial/binarysplitfact.html">http://www.luschny.de/math/factorial/binarysplitfact.html</a>"<br>    p = 1<br>    r = 1<br>    p, r = factorial_loop(n, p, r)<br>
    return r << n_minus_num_of_bits(n)<br><br>def factorial_loop(n, p, r):<br>    if n <= 2: return p, r<br>    p, r = factorial_loop(n // 2, p, r)<br>    p *= part_product(n // 2 + 1 + ((n // 2) & 1),<br>                      n - 1 + (n & 1))<br>
    r *= p<br>    return p, r<br><br>def part_product(n, m):<br>    if (m <= (n+1)): return n<br>    if (m == (n+2)): return n*m<br>    k = (n+m) // 2<br>    if ((k & 1) != 1):<br>        k -= 1<br>    return part_product(n, k) * part_product(k+2, m)<br>
<br>def n_minus_num_of_bits(v):<br>    w = v<br>    if w >= 2**32-1:<br>        raise OverflowError<br>    w -= (0xaaaaaaaa & w) >> 1<br>    w = (w & 0x33333333) + ((w >> 2) & 0x33333333)<br>    w = w + (w >> 4) & 0x0f0f0f0f<br>
    w += w >> 8<br>    w += w >> 16<br>    return v - (w & 0xff)<br></div></div><blockquote style="margin: 1.5em 0pt;">--<br>
Daniel Stutzbach, Ph.D.<br>
President, <a href="http://stutzbachenterprises.com">Stutzbach Enterprises, LLC</a>
</blockquote>