Beefy Boxes and Bandwidth Generously Provided by pair Networks
Pathologically Eclectic Rubbish Lister
 
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??


in reply to Re^2: Fastest way to calculate hypergeometric distribution probabilities (i.e. BIG factorials)?
in thread Fastest way to calculate hypergeometric distribution probabilities (i.e. BIG factorials)?

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)
  • 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 ]; } }

    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) ^

    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 browsing the Monastery: (10)
As of 2019-11-18 16:37 GMT
Sections?
Information?
Find Nodes?
Leftovers?
    Voting Booth?
    Strict and warnings: which comes first?



    Results (90 votes). Check out past polls.

    Notices?