Clear questions and runnable codeget the best and fastest answer PerlMonks

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

by salva (Abbot)
 on Jun 14, 2005 at 14:02 UTC ( #466542=note: print w/replies, xml ) Need Help??

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)
[download]```
• Comment on Re^3: Fastest way to calculate hypergeometric distribution probabilities (i.e. BIG factorials)?
• Download Code

Replies are listed 'Best First'.
Re^4: Fastest way to calculate hypergeometric distribution probabilities (i.e. BIG factorials)?
by tlm (Prior) on Jun 16, 2005 at 00:56 UTC

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 ];
}
}
[download]```

the lowliest monk

Re^4: Fastest way to calculate hypergeometric distribution probabilities (i.e. BIG factorials)?
by BrowserUk (Pope) on Jun 14, 2005 at 14:35 UTC

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

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

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

Examine what is said, not who speaks -- Silence betokens consent -- Love the truth but pardon error.
Lingua non convalesco, consenesco et abolesco. -- Rule 1 has a caveat! -- Who broke the cabal?
"Science is about questioning the status quo. Questioning authority".
The "good enough" maybe good enough for the now, and perfection maybe unobtainable, but that should not preclude us from striving for perfection, when time, circumstance or desire allow.

Log In?
 Username: Password:

What's my password?
Create A New User
Node Status?
node history
Node Type: note [id://466542]
help
Chatterbox?
and the web crawler heard nothing...

How do I use this? | Other CB clients
Other Users?
Others examining the Monastery: (3)
As of 2020-10-28 09:53 GMT
Sections?
Information?
Find Nodes?
Leftovers?
Voting Booth?
My favourite web site is:

Results (260 votes). Check out past polls.

Notices?