P is for Practical 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));
}

+
```PURE     : 40320
STIRLING : 417351.253768542
REAL     : 39902.3954526567
GOSPER   : 40034.4823215513

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 browsing the Monastery: (4)
As of 2020-01-27 10:32 GMT
Sections?
Information?
Find Nodes?
Leftovers?
Voting Booth?
The worst excuse I have ever heard is:

Results (112 votes). Check out past polls.

Notices?