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.

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.

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';
[download]

I can use this with the normal computation

use Math::GMP qw(:constant);
[download]

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.


Regards,
Edward

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

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

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%       --
[download]

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

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

the lowliest monk

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));
[download]

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

That formula has an error. It should be

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

log(exp(-n)) = -n

-Mark

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

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:

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

thor

Feel the white light, the light within
Be your own disciple, fan the sparks of will
For all of us waiting, your kingdom will come

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