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

[jonadab]

though I don't fully understand at what size number these modules become necessary

*Much* smaller numbers than the ones you're dealing with. The details are platform-dependent; for instance, a 64-bit platform can handle, without a big-number library, larger numbers than a 32-bit platform. But nobody is ever going to manufacture a computer that can natively handle numbers of the size you're talking about. Based on Moore's law, you couldn't expect such a computer to be available while Earth is still inhabitable.

You may want to look into ways to approximate factorials. Irrespective of what number library you use, I'm not sure it's possible to calculate the exact value of the factorial of 700 in reasonable time on a household computer.

Alternatively, look for ways to shortcut the problem, by simplifying before you multiply. You may find that some things cancel out in ways that save a lot of time. Also, if this is for a math class, you may want to check, as some math teachers are quite happy to receive answers in expression form, with not all of the calculations performed, because they're more interested in whether you understand how to work the problem than in whether you can multiply together a few hundred numbers.

[BrowserUk]

What do you consider "a reasonable time"?

P:\test>perl -mTime::HiRes=time -MMath::BigInt 
-wle"print time; $n=Math::BigInt->new( 700 ); print $n->bfac; print ti
+me;"
1118754720.06049
2422040124750272179867875093812352218590983385729207299450679664929938
1602156474204445190516664848192493214566714970498423275250938748173438
3839375763145922845082849997227127414016031105783055846363633712407933
2447820739281101037112665387537180790257577919273108262916904750405235
0550600840122194928923756351362966220200231781096198180461799068974504
2054891261087058908805650391358456221103769328878296090019507413099979
9035970711436279339094292032866260496375825461427727555710003007752906
1414706395743900249885149142644498650064588732269519418995459703339103
5158855923294082956927698608022220028996612834393163002878920338265474
9603473516314765262772257171154686716862814184728741187147936349501653
1974574556604131345060491220449470526233846820888647906733095692923842
1561178801427495490591414836230322620024681644130193484608025499864732
5270606104512088058712293349862185399243309054299576381718806247238195
2326046426143298940706361637536720912327516123783482738407578735677175
3287924251833711954060294360941162934900956604372083673740109088239297
5031224612531245642687296717053747734506443314924558119560479901478736
2095569251615177371103997547305518540663284200147286578962869365237870
8020647632715713644131877343275100726310805695825169381128095724320246
0157111778617472683761623869704457588005158037495665069625778930898095
7257947107016392382315281155796191202873786892389343351985086659339172
5714397527770759059751198934506870173594016967256186471310711501674736
8992690116082633762172346688969840862517264384000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000
1118754720.13248

[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.