Beefy Boxes and Bandwidth Generously Provided by pair Networks
laziness, impatience, and hubris
 
PerlMonks  

Re: Stirling Approx to N! for large Number?

by tall_man (Parson)
on Mar 25, 2005 at 15:46 UTC ( #442358=note: print w/ replies, xml ) Need Help??


in reply to Stirling Approx to N! for large Number?

There was a more accurate formula than Stirling's in tilly's post. Wolfram's Mathematica attributes it to Gosper:

sub factorial_gosper { # sqrt(2*n*PI + 1/3)*n**n/exp(n) my $n = shift; return exp($n*log($n) - $n + 0.5*log(2*$n*PI + 1.0/3.0)); }

This gives better accuracy for small numbers. For 8!, I get 40034.48, compared to the correct 40320.

There are several problems with your implemenation of Stirling's formula. For one thing, log(e^(-n)), which you give as -e*log(n). Since the log is just the exponent base e, the correct answer is -n.

Tilly's formula is closer to being right; there was only one parenthesis mistake. I corrected that and tried them all in a modified version of your code:

#!/usr/bin/perl -w use strict; + my $h =8; my $pure = factorial_pure($h); my $stirling = factorial_stirling($h); my $real_stirling = real_stirling($h); my $gosper = factorial_gosper($h); + print "PURE : $pure\n"; print "STIRLING : $stirling\n"; print "REAL : $real_stirling\n"; print "GOSPER : $gosper\n"; + #------Subroutines------------ + sub factorial_stirling{ my $n = shift; use constant PI => 4*atan2 (1,1); use constant e => exp(1); + my $log_nfact = $n * log ($n) - e* log($n) + 0.5 * (log(2*PI)+log($n)); + # Below is Tilly's suggestion, with less accuracy # $n * log ($n) - ($n) + (0.5 * (log(2+ log(PI))+log($n))); + return exp($log_nfact); } + + sub factorial_pure { + my ($n,$res) = (shift,1); return undef unless $n>=0 and $n == int($n); $res *= $n-- while $n>1; return $res; + } + sub real_stirling { my $n = shift; my $log_nfact = $n*log($n) - $n + 0.5*(log(2) + log(PI)+ log($n)); return exp($log_nfact); } + sub factorial_gosper { # sqrt(2*n*PI + 1/3)*n**n/exp(n) my $n = shift; return exp($n*log($n) - $n + 0.5*log(2*$n*PI + 1.0/3.0)); } +
I get the answers below:
PURE : 40320 STIRLING : 417351.253768542 REAL : 39902.3954526567 GOSPER : 40034.4823215513


Comment on Re: Stirling Approx to N! for large Number?
Select or Download Code

Log In?
Username:
Password:

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

How do I use this? | Other CB clients
Other Users?
Others making s'mores by the fire in the courtyard of the Monastery: (17)
As of 2015-07-30 20:39 GMT
Sections?
Information?
Find Nodes?
Leftovers?
    Voting Booth?

    The top three priorities of my open tasks are (in descending order of likelihood to be worked on) ...









    Results (273 votes), past polls