Originally posted by CP5670
Do you mean Stirling's formula? (Binet's formula is the bernoulli-type e2pt-1 integral) Although the Lanczos approximation is also quite effective, since it is derived from the same thing except with a term dropped out. 
There is a pretty nice set of numerical special functions out on the internet for the 89 and 92+. Includes gamma, zeta, legendre P and a few others. It cannot do symbolic work, but is good for getting graphs and stuff of the functions.
I'm pretty sure it was a "form" of the Lanczos.
Reference:
http://www.rskey.org/gamma.htmLook around where it says "...demonstrates a much more efficient algorithm, the so-called Lanczos approximation for computing the Gamma function..." and then you'll start running into:
Gamma(z) = Sqrt(2*pi)/z*(P(0) + [n = 1 through 6] P(n)/(z+n)) * (z+5.5)^(z+0.5)*e^-(z+5.5)
where:
p0 = 1.000000000190015
p1 = 76.18009172947146
p2 = -86.50532032941677
p3 = 24.01409824083091
p4 = -1.231739572450155
p5 = 1.208650973866179 * 10^(-3)
p6 = -5.395239384953 * 10^(-6)
And then from that:
Gamma(-z) = -pi/(Gamma(z+1)*Sin(pi*z))
And using the recurrence relation, the factorial of much larger numbers can be computed (so far for one of my java programs I got the Gamma(171) = 7.257415615308e306 where Windows Calculator computed 170! to be 7.257415615307998967396728211129e+306 so I'd say this thing is pretty accurate
).
If you want my source code and stuff just gimme a hollar.
There's lots of other good stuff on that page also.
Edit: Might I note that the Stirling formula wasn't "precise" enough for what I was trying to calculate (the shorter one, anyway as shown on the page... and for some reason the long version didn't quite help me out either, not sure why
). Although x^x*e^-x*Sqrt(2*pi*x)*(1+1/(12*x)) looks pretty good on a graph, the number output for the approximation above looked much, much better.
Edit #2: Oh yea, I remember why the longer version of Stirling's wasn't working right - it was because the accuracy only improved for digits over 5.
This one?:
e^(x*ln(x)-x+ln(Sqrt(2*3.14159265358979/x))+1/(1188*x^9)-1/(1680*x^7)+1/(1260*x^5)-1/(360*x^3)+1/(12*x))
