This is more than an order of magnitude faster which should get you into the realms of web-ability.

Your code takes 0.000063938 seconds; the code below 0.000002392 seconds; which should reduce your 70 seconds to 2.6.

I've precalculated all the logfacts for 0 .. 1000. If your inputs need a different range, recalculate:

#!/usr/bin/perl -slw
use strict;

my @logfact = <DATA>;

sub hypergeom {
    # There are m "bad" and n "good" balls in an urn.
    # Pick N of them. The probability of i or more successful selectio
+ns:
    # (m!n!N!(m+n-N)!)/(i!(n-i)!(m+i-N)!(N-i)!(m+n)!)
    my ($n, $m, $N, $i) = @_;

    my $loghyp1 = $logfact[ $m ] + $logfact[ $n ]
                + $logfact[ $N ] + $logfact[ $m + $n - $N ];
    my $loghyp2 = $logfact[ $i ] + $logfact[ $n - $i ] 
                + $logfact[ $m + $i - $N ] + $logfact[ $N - $i ] 
                + $logfact[ $m + $n ];

    return exp($loghyp1 - $loghyp2);
}

print hypergeom( 300, 700, 100, 40 );
[download]

sub hypergeom1 {
    # There are m "bad" and n "good" balls in an urn.
    # Pick N of them. The probability of i or more successful selectio
+ns:
    # (m!n!N!(m+n-N)!)/(i!(n-i)!(m+i-N)!(N-i)!(m+n)!)
    my ($n, $m, $N, $i) = @_;

    my $loghyp1 = logfact1($m) +logfact1($n)+logfact1($N)+logfact1($m+
+$n-$N);
    my $loghyp2 = logfact1($i)+logfact1($n-$i)+logfact1($m+$i-$N)+logf
+act1($N-$i)+logfact1($m+$n);
    return exp($loghyp1 - $loghyp2);
}

sub logfact1 {
    my $x = shift;
    my $ser = (   1.000000000190015
                + 76.18009172947146   / ($x + 2)
                - 86.50532032941677   / ($x + 3)
                + 24.01409824083091   / ($x + 4)
                - 1.231739572450155   / ($x + 5)
                + 0.12086509738661e-2 / ($x + 6)
                - 0.5395239384953e-5  / ($x + 7) );
    my $tmp = $x + 6.5;
    ($x + 1.5) * log($tmp) - $tmp + log(2.5066282746310005 * $ser / ($
+x+1));
}
[download]

#!/usr/bin/perl 

use strict;
use warnings;

sub logfact {
   return gammln(shift(@_) + 1);
}

sub hypergeom {
   my ($n, $m, $N, $i) = @_;
   my $loghyp1 = logfact($m)+logfact($n)+logfact($N)+logfact($m+$n-$N)
+;
   my $loghyp2 = logfact($i)+logfact($n-$i)+logfact($m+$i-$N)+logfact(
+$N-$i)+logfact($m+$n);
   return exp($loghyp1 - $loghyp2);
}

sub gammln {
  my $xx = shift;
  my @cof = (76.180091729471457, -86.505320329416776,
             24.014098240830911, -1.231739572450155,
             0.12086509738661, -5.395239384953e-06);
  my $y = my $x = $xx;
  my $tmp = $x + 5.5;
  $tmp -= ($x + 0.5) * log($tmp);
  my $ser = 1.0000000001900149;
  foreach my $j (0 .. 5) {
     $ser += $cof[$j] / ++$y;
  }
  -$tmp + log(2.5066282746310007*$ser/$x);
}

print hypergeom(300,700,100,40),"\n";
[download]

Minor nit: As the author of Math::GammaFunction, I would consider Math::GSL a more mature, more reliable implementation.

--Steffen

mrguy123:

Have you tried to use Memoize to cache the results of logfact? If it's repeatedly being called with the same values, you might get a speed boost. You could try to Memoize hypergeom as well--but if there are too many combinations of arguments (i.e., the argument list rarely repeats) it would be a loss.

...roboticus

Update: Just for my own amusement, I added Memoize to the program and found no significant improvements. Apparently the functions are fast enough that there's no big win to be found simply by caching the values. I replaced the single call at the end of the program with a nested loop (to generate enough runtime to be significant). This should call hypergeom 203,401 times:

for my $i (280 .. 320) {
    for my $j (680 .. 720) {
        for my $k (95 .. 105) {
            for my $l (30 .. 40) {
                my $hg = hypergeom($i, $j, $k, $l);
            }
        }
    }
}
[download]

The results are:

39.11s   Without Memoize
35.25s   Memoize logfact only
38.40s   Memoize gammaln only
47.04s   Memoize hypergeom only
[download]

Of course, with the usage pattern given, I may have artificially improved/reduced Memoize's results, so your mileage will vary.

PDL is very fast, as fast as C when it is number crunching. Google for "PDL hypergeometric" for ideas. The PDL maillist is very helpful.

I have an input set that needs to call the hypergeom function 137544 times

And just as an aside: if installing modules locally is a problem due to not having root access: check out local::lib, which allows you to do just that.