Beefy Boxes and Bandwidth Generously Provided by pair Networks
Keep It Simple, Stupid
 
PerlMonks  

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

by BrowserUk (Pope)
on Jun 14, 2005 at 13:15 UTC ( #466517=note: print w/replies, xml ) Need Help??


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

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

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

Replies are listed 'Best First'.
Re^3: Fastest way to calculate hypergeometric distribution probabilities (i.e. BIG factorials)?
by salva (Abbot) on Jun 14, 2005 at 14:02 UTC
    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)

      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

      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://466517]
help
Chatterbox?
and the web crawler heard nothing...

How do I use this? | Other CB clients
Other Users?
Others pondering the Monastery: (7)
As of 2020-01-28 19:11 GMT
Sections?
Information?
Find Nodes?
Leftovers?
    Voting Booth?
    Notices?