On Sat, Jan 23, 2010 at 4:55 PM, kj <no.em...@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>
-- http://mail.python.org/mailman/listinfo/python-list
