Stirling Approx to N! for large Number?

[monkfan]

Hi,
Following tilly's posting here. I attempted to construct the code.

#!/usr/bin/perl -w
use strict;

my $h =8;
my $pure     = factorial_pure($h);
my $stirling = factorial_stirling($h);

print "PURE     : $pure\n";
print "STIRLING : $stirling\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;

    }

Stirling method indeed gives a faster result:

             Rate     pure stirling
pure     290034/s       --     -30%
stirling 319835/s      16%       --

With a reasonable approximation:

PURE     : 40320
STIRLING : 417351.253768542

But when it comes to large N, e.g $h = 500. It began to fail.

I wonder if Tilly's suggestion is ever be possible,
i.e resorting to pure Perl arithmetic (no external module) yet able to handle large number in a fastway?

[Zaxo]

What's pure? Perl core has Math::BigInt which can use GMP or Pari to speed up representation and operations.

Pure usually means "written in Perl" but linkage to libraries is not immoral or dishonest.

use Math::BigInt qw(lib GMP :constant);

Update: corrected :constant tag. (Later, corrected with 'lib' in response to ewijaya's reply). The :constant tag is not wanted if you are dealing with floats.

[monkfan]

Zaxo,
What's the equivalent to that in Pari? I tried this at the beginning of my code. But wont do.

use Math::BigInt lib => 'Pari';

I can use this with the normal computation

use Math::GMP qw(:constant);

But when its applied with my "factorial_stirling" subroutine, it became in conflict with my 'use constant statement for PI, EXP and also LOG function.

[Fletch]

I hope you miscopied your "reasonable approximation" because it's off by an order of magnitude . . . (4.0e5 vs 4.1e6). Also IMSMR log(n*m) should be log(n) + log(m) so the commented out example should be log(2) + log(PI) rather than log(2 + log(PI)) (which might explain any inacurracy).

Addendum: If you have access to a copy of Numerical Recipies In C (ISBN 0521431085), chapter 6 has a factorial approximation using ln(Γ(x)) which should be easily translated into Perl (or use it unaltered via Inline::C).

[jimbojones]

Hi
Numerical Recipes in C is available free online. See Numerical Recipes Books On-Line. It's a great resource. The site says:

Thanks to special permission from Cambridge University Press, we are able to bring you, free, the complete Numerical Recipes books in C, Fortran 77, and Fortran 90 On-Line, in Adobe Acrobat format. Due to copyright restrictions, Numerical Recipes in C++ is not available as part of this free resource.

The factorial representation is available at this PDF (after a I chose LANL as the mirror site).
- j

[tlm]

Factorial is a good candidate for memoization. A memoized version of factorial is faster than Stirling's (at least for a naive implementation of Stirling's in the range of arguments that don't cause overflow). Here are the stats:

             Rate     pure stirling     memo
pure      15037/s       --     -89%     -94%
stirling 139183/s     826%       --     -48%
memo     267957/s    1682%      93%       --

I give the code below. Note that I fixed your Stirling approximation; I used one of the versions given by MathWorld (expression 11). I also made some other small changes.

On my system, all these subroutines overflow for any argument > 170.

Afterthought: Actually, this range is so small that you may as well save yourself the subroutine call: just precompute the whole array, and access $fact[ $n ] directly. On the other hand, retaining the procedural interface is better software engineering.

Update 2: Fixed the handling of invalid arguments in factorial_memo.

the lowliest monk

Read more... (2 kB)

[tilly]

With the way that Perl chooses to handle numbers, you cannot represent large numbers in native Perl. (Other languages do not necessarily have this limitation.) Anyways a correction and a suggestion.

The correction is that you misunderstood what I suggested would be more accurate. Try this line for the log(n!):

my $log_nfact =
  $n * log ($n) - $n + 0.5 * (log(2*PI*$n + 1/3));

Try it, it should be better.

The suggestion is that your problem was about binomial probabilities. Rather than try to calculate the factorials and work with them directly, estimate the logs of the factorials and use them to calculate the log of n choose m. Only then use the exponential.

For n and m large enough you'll again run into problems with Perl's support for large numbers. But since the final number is more modest than the intermediate factorials, you'll be able to go much farther. And when you get beyond what Perl can handle natively, if you have the log, you can divide it by log(10) to get the log base 10, and now you can take the integer and fractional parts of that and readily convert your answer into scientific notation for display even though Perl's native scientific notation doesn't go that high.

And that latter strategy will work for any n and m that you're likely to encounter in practice.

Update: Fixed algebra per kvale's note. (What I get for posting fast on a late night.)

[kvale]

That formula has an error. It should be

my $log_nfact =
  $n * (log ($n) - 1) + 0.5 * (log(2*PI*$n + 1/3));

log(exp(-n)) = -n

[tall_man]

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:

Read more... (4 kB)

[chas]

"$n * log ($n) - e* log($n) + 0.5 * (log(2*PI)+log($n));"
Stirling's approx is n^ne^{-n}sqrt{2 pi n), so when you take the log, you get nlog(n) - nlog(e) + .5(log(2pi)+log(n)) so it appears one of your terms is incorrect. (Why not use natural log anyway?)
In addition Stirling's gives an asymptotic approximation (the ratio of Stirling's expression to n! is close to 1 for large n; it is not true that the expresion is close to n! for large n.)
chas

[thor]

Why not use natural log anyway?

In Perl, the log function does return the natural log. However, it is a property of logarithms that the logarithms of the same number in a different base differ only by a constant multiplicitave constant. That is to say that log_a(x) = log_b(x)/log_b(a). So if one person used their favorite base and you have a different favorite base, you need only scale their answer by a factor of log_yours(theirs).

[chas]

Yes, I guess I was thinking the OP meant $n*log(e) rather than e*log($n) (the former is correct), but then log(e) = 1 if log is the natural log (so I wondered if some other log was meant.) You're correct that my comment about using natural log is extraneous...thanks.
chas