Re^3: Fastest way to calculate hypergeometric distribution probabilities (i.e. BIG factorials)?

[salva]

or when absolute precission is not required, and in pure perl:

#!/usr/bin/perl

use Time::HiRes 'time';

my $next = 1;
my @log;
my @fact=(0);

sub more_logs {
    my $n = shift;
    for ($next..$n) {
    my $log =log($_);
    $log[$_] = $log;
    $fact[$_] = $fact[$_-1]+$log;
    }
    $next = $n+1;
}

more_logs(10);
my $l10 = $log[10];
my $il10 = 1.0/$l10;

for (@ARGV) {
    my $start = time;
    more_logs($_) if $_ >= $next;
    my $fact = $fact[$_];
    my $exp = int $fact*$il10;
    my $man = exp($fact-$exp*$l10);
    my $end = time;
    printf("%d! = %0.10fE%d (time: %f-%f)\n", $_, $man, $exp, $start, 
+$end)
}

__END__ 

$ perl /tmp/bigfact.pl 1000 700 600 2000 5000
1000! = 4.0238726008E2567 (time: 1118757636.805467-1118757636.810242)
700! = 2.4220401248E1689 (time: 1118757636.810742-1118757636.810792)
600! = 1.2655723162E1408 (time: 1118757636.810928-1118757636.810945)
2000! = 3.3162750925E5735 (time: 1118757636.811015-1118757636.816799)
5000! = 4.2285779269E16325 (time: 1118757636.816961-1118757637.117728)

[tlm]

Nice, but neither @log nor $next are necessary, and the whole thing can be given a functional interface, while transparently memoizing the results:

{
  my @ln_fact;
  BEGIN { $ln_fact[ 0 ] = 0; }

  sub ln_fact {
    # returns the natural log of the factorial

    my $n = shift;
    die "Invalid arg: $n" if $n < 0 or $n > int $n;
    $ln_fact[$_] = $ln_fact[$_-1] + log($_) for @ln_fact..$n;
    return $ln_fact[ $n ];
  }
}

[BrowserUk]

S'neat, but you do loose quite a bit more accuracy than with Math::Pari

Math::BigInt 19937! = 2.174930234150431374566653562e77066 time: 104.69
[phone]
Math::Pari   19937! = 2.174930234150431374566653166e77066 time:   0.00
[phone]
                                                ^
PP::LogFact  19937! = 2.174930232142160900000000000e77066 (time:  0.03
+0866)
                                ^

Once you start doing further math with them, those inaccuracies add up fast.