From: casevh on 24 Jan 2010 16:49 On Jan 23, 2:55 pm, kj 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. > > TIA! > > ~K If you need exact values, gmpy (http://code.google.com/p/gmpy/) has basic, but fast, support for calculating binomial coefficients and factorials. If you are floating point approximations, check out mpmath (http://code.google.com/p/mpmath/). casevh From: Mel on 24 Jan 2010 19:03 Dave Angel wrote: > But I still think there must be code for the gamma function that would > be quicker. But I haven't chased that lead. Log of the gamma function is also an algorithm in its own right (e.g. published in _Numerical Recipes_.) Mel. From: Peter Pearson on 27 Jan 2010 23:33 On Sun, 24 Jan 2010 16:11:03 -0500, Dave Angel wrote: > I didn't think of simply summing the logs. A couple terms of Stirling's approximation work pretty well: def log_fact_half( N ): """log_fact_half( n ) returns the natural logarithm of the factorial of n/2. n need not be an integer. Domain: n from 0 to infinity Range: from something around -0.12 to infinity """ SmallFacts = ( 0.0000000000000000, -0.1207822376352453, 0.0000000000000000, 0.2846828704729192, 0.6931471805599453, 1.2009736023470743, 1.7917594692280550, 2.4537365708424423, 3.1780538303479458, 3.9578139676187165, 4.7874917427820458, 5.6625620598571418, 6.5792512120101012, 7.5343642367587327, 8.5251613610654147, 9.5492672573009969, 10.6046029027452509, 11.6893334207972686, 12.8018274800814691, 13.9406252194037634 ) if N < 0: raise RuntimeError, \ "Negative argument in LogHalfFact!" if N < len( SmallFacts ): RetVal = SmallFacts[ N ] else: n = 0.5 * N RetVal = (n+0.5)*math.log(n) - n + HALFLOG2PI + 1.0/(12*n) - \ 1.0/(360*n*n*n) return RetVal -- To email me, substitute nowhere->spamcop, invalid->net.