ISO module for binomial coefficients, etc.

[kj]

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

 

[idrevetnom]

Maybe this could be of interest :


http://tnt.math.metro-u.ac.jp/nzmath/manual/modules/combinatorial.html


hope this helps Id

 

[Dave Angel]

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

You do realize that a standard. python floating point number cannot
possibly approximate a number like 10000! Better use longs.

I'd check out the gamma function, which matches factorial for integer
arguments (plus or minus 1).

DaveA

 

[Alf P. Steinbach]

* Dave Angel:

You do realize that a standard. python floating point number cannot
possibly approximate a number like 10000!

I think what kj is looking for, provided she/he is knowledgable about the
subject, is code that does something like
... log_fac += log( i, 10 )
...
which turned out to be accurate to 10 digits.

Better use longs.

That would involve incredible overhead. E.g., how many bytes for the number
above? Those bytes translate into arithmetic overhead.

I'd check out the gamma function, which matches factorial for integer
arguments (plus or minus 1).

Or, e.g., logarithms... ;-)


Cheers & hth.,

- Alf

 

[Dave Angel]

I think what kj is looking for, provided she/he is knowledgable about
the subject, is code that does something like

... log_fac += log( i, 10 )
...
) ) )
10000! = 2.84625968062e35659

which turned out to be accurate to 10 digits.



That would involve incredible overhead. E.g., how many bytes for the
number above? Those bytes translate into arithmetic overhead.
About 14k.


Or, e.g., logarithms... ;-)


Cheers & hth.,

- Alf

I didn't think of simply summing the logs. I did have some
optimizations in mind for the multiply of the longs. If you do lots of
partial products, you can do a good bit of the work with smaller
numbers, and only get to longs when those partial products get big
enough. You could also use scaling when the numbers do start getting
bigger.

But I still think there must be code for the gamma function that would
be quicker. But I haven't chased that lead.

DaveA

 

[casevh]

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

 

[Mel]

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.

 

[Peter Pearson]

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