Also isn't exponentiating the log of n! and carrying accuracy to the last integer is as much or more work than calculating n! in the first place? If you just want the first k digits, the approximation is better. Most of us don't.

That said, the approximations are good for order of magnitude estimates. Understanding the basic approximations is easier than most think. Watch:

log(n!) = log(1) + log(2) + ... + log(n)
        = 1/2 log(1) + (log(1) + log(2))/2
                     + (log(2) + log(3))/2
                     + ...
                     + (log(n-1) + log(n))/2
                     + 1/2 log(n)
[download]

log(n!) = 1/2 0 + Integral from 1 to n of log(x) +
          1/2 log(n) + e(1) + e(2) + ... + e(n-1)
[download]

The answer is, "Quite a bit." The goal of approximation techniques is to take something, divide it into a few big terms, understand those terms, a big mess, show that the mess isn't important. This is what we have done.

That integral is a single huge term. But it is one we understand. The integral of log(x) is x*log(x)-x. So that integral is n*log(n)-n-(1*log(1)-1) = n*log(n)-n+1. log(n)/2 is the next biggest thing around. And finally we have those e()'s. We think that they are small. But how to show that? It turns out that the error in a trapezoidal approximation is w**3 * f''(t)/12 where w is the width of the interval and t is a point in it. (This is reasonable - the error in a linear approximation is tied to what the function's bend does to the integral.) The second derivative of log(x) is -1/(x*x), so e(i) is O(1/(i*i)). If you remember your infinite series, that means that the infinite sum e(1)+e(2)+e(3)+... converges to some number, call it E. And the partial sum from n out to infinity is O(1/n) from that number.

Remember the principle of approximation? Get a few big understood terms, at the cost adding lots of messy small ones? We are about to do that where "lots" is an infinite number! Taking the above we get:

log(n!) = 1/2 0 + Integral from 1 to n of log(x) +
          1/2 log(n) + e(1) + e(2) + ... + e(n-1)
        = n*log(n) - n + 1 - 1/2 log(n)
          + (e(1) + e(2) + e(3) + ...)
          - (e(n) + e(n+1) + e(n+2) + ...)
        = n*log(n) - n + log(n)/2 + E + 1
          - e(n) - e(n+1) - e(n+2) - ...
[download]

n! = K * sqrt(n) + n**n * (1 + O(1/n)) / e**n
[download]

The first is straightforward but messy. Just look at the infinite sum as an approximation to an integral. Write it as that integral and associated error terms. Get another term of the approximation, with smaller and more complex error terms. Wash rinse and repeat to your heart's desire.

The second requires cleverness. I suggest sticking this approximation into the binomial expansion, and compare with the central limit theorem's approximation to show that K has to be 1/sqrt(2*pi). If you are clever enough about it you can actually derive the central limit theorem for the binomial expansion up to the constant, then show that the only value of K which makes the terms for (1+1)**n come out to 2**n for large n is 1/sqrt(2*pi).

These are left as exercises because I am tired of typing. :-)

Why do you think that truncating Stirling's series is preferable to approximating the more accurate Lanczos asymptotic series (which is valid for Re(n) > 0)? It is more 'robust' to calculate the logarithm of the gamma function because you are less prone to overflow for large values of n.

I didn't say it was preferable. I tried to explain how you can derive one of these approximations.

As for robustness, all of the solutions presented in this discussion use big integer arithmetic. Until you run out of memory, they are precise to the last digit.