Prev: how to generate random numbers that satisfy certain distribution
Next: image processing - inverse filtering
From: casevh on 24 Jan 2010 16:49
On Jan 23, 2:55 pm, kj <no.em...(a)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.
>
> 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 <davea(a)ieee.org> 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.
First  |  Prev  | 
Pages: 1 2
Prev: how to generate random numbers that satisfy certain distribution
Next: image processing - inverse filtering