An Asymptotic Derivation
The above nodes only give a single term - but we can get an infinite series if we try slightly harder. This series will be asymptotic and diverge - but will still be very useful if we only take the first few terms. Let us consider what we actually mean by a factorial. Well, it can be shown that the factorial function (defined on the natural numbers) has precisely one analytic continuation given suitable assumptions about convexity or equivalent - the Gamma function.
Gamma (z) == ∫ exp (-t + (z-1)ln (t)) dx
Where the integral runs between 0 and infinity along the real line. Take a look at that integrand. It's largest when t = (z-1), and is pretty damned small everywhere else. I'm going to apply Laplace's Method to the region around t = (z-1) and see if I can't get somewhere1.
Write h(t) === -t + (z-1)ln (t) and Taylor expand around t=(z-1).
h(t) ~ (z-1)( -1 + log (z-1)) + 0*t - (t - (z-1))*(t - (z-1))/(2*(z-1)) + O((t - (z-1))^{3})
Why have I stopped there? After all, t is a variable - it runs from 0 to infinity, so I can hardly claim it's small. Well, it's not small. But h(t) is quite negative away from t = (z-1), and so the integrand is exponentially small; so, for my leading order term, I'm going to stop here. If I wanted higher order terms, I would simply keep more terms from h.
Okay, so this gives us
Gamma(z) ~ (z-1)^(z-1) * exp(1-z) * ∫ exp( - u*u) du*(sqrt(2 * (z - 1))
Where I have sneakily made the substitution (t-(z-1))*(t-(z-1))/(2*(z-1)) = u*u, and the limits of the integral are now from -infinity to + infinity. Remember that away from t = (z-1), everything is pretty much zero (those magic words were exponentially small) and it doesn't matter if we integrate over that space or not. It's convenient if we do. We recall our favourite Gaussian integral is equal to sqrt(pi) and that, for this to be a good approximation, z must be large, to find:
z! ~ Gamma(z+1) ~ z^(z + 1/2) * sqrt(pi) exp(-z).
Higher Order Terms
There are two approximations I made: To neglect the exponentially small region (which is really, really, really tiny, okay? There's no point keeping it here - although that's not always true) and to drop some of the larger terms from h. If I'd kept them in h, then I would have ended up with something that looks like:
exp(A + Bu*u + Cu*u*u + ...)
in my integrand. Well, an easy thing to do here (though by no means the only thing) is to expand out some of that exponential:
=exp(A + B u*u + 0 + ...)(1 + Cu*u*u + ...)
which is another integral I can evaluate, term by term. Note that with the substitution I made earlier I'd pick up another power of z^{-1/2}, making this term in some sense smaller than the previous. Then I'd find a term which looked like 1/(12*z) in my expansion; the next lowest term.
Why do this?
This one? It's a nice example of Laplace's Method, and is one of the more basic asymptotic problems. In general, if one has a nasty integral which one cannot evaluate analytically, it is often convenient to find an asymptotic expansion which is cheap computationally and retains the basic mathematical ingredients. Please believe me when I say this is not some cheap trick one might use once, but is a powerful tool in applied mathematics.
See Also