An approximation to the factorial function.

N! ~ N^N e^-N

or, in logarithmic form:

ln N! = N ln N - N

For large N, this is a good approximation to the factorial.

With some adjustments it is possible to make it better, but since it is mostly used for back-of-the-envelope types of calculations, this is usually good enough.

This is based on a truncation of Stirling's Series, an infinite series which converges to exactly the value of the factorial.

From the standpoint of a number theorist, Stirling's formula is a significantly inaccurate estimate of the factorial function (n! = 1*2*3*...*(n-1)*(n)). The formula is:

n! ~ (n/e)ⁿ

There are a couple ways of deriving this result. One is to use the Euler Gamma function; take the derivative of the functional portion of the integral, then work from there. Another is to use the natural logarithm directly:

ln(n!) = ln (n*(n-1)*(n-2)*...*2*1) = ln n + ln (n-1) + ... + ln 2 + ln 1.

This is not an approximation, but it is also functionally equivalent to computing n! directly, so it is equally useless. However:

           n           / n
         -----        /
         \           |
ln(n!) =  > ln(i) ~   \   ln(x)dx = nln(n)-n+1.
         /             |
         -----        /
          i=1        / 1

e^ln(n!)=n!, so n! ~ e^nln(n)-n+1=e*(n/e)ⁿ. We leave off the 1 because, as n approaches infinity, the 1 "becomes insignificant". For small values of n, the formula is wildly inaccurate regardless of whether the 1 is present.

Why on earth would I want to use this? n! can be calculated directly, right?

True. For small values of n, the best way to get n! is to calculate it directly. But what does it mean for n to be "small"? That probably depends on the computing equipment at hand. Is it particularly time intensive to calculate 300!? 3000!? 3000000!? How about 10⁴⁰!?

As of the timestamp of this writing, the largest known prime number is the Mersenne prime 2^13,466,917-1, which has roughly 4 million decimal digits. From Wilson's theorem, the fastest way to check if a number is prime might be to use the congruence -1=(n-1)! (mod n) (since there is only one test involved). So, can your TI-89 calculate the value of (2^13,466,917-2)!? I wish it could. Because of the sheer size of the numbers involved, an approximation of n! would be very nice to have, if only it were more accurate. Another significant block in the usage of this formula is the computational accuracy of the number e. nⁿ is relatively easy, but eⁿ=1+n+(n²)/2!+(n³)/3!+... is not accurate at the integer level before the nth term (nⁿ)/n! (nor for quite a distance after), which involves the value in question. Thus, as it stands, Stirling's formula is merely of passing interest in number theory.

---

Note: the formula that krimson states is the "real" formula that Stirling proposed; it's an approximation given by truncating all but the first term of the Stirling Series (see http://mathworld.wolfram.com/StirlingsApproximation.html, where other even closer approximations are given). The point of my writing here is to make people aware of how useful the formula is in other contexts.

[The "significant inaccuracy" is dramatically reduced by using the correct formula]

n! ~ (2πn)^1/2nⁿe^-n

Stirling's formula provides an approximation to the factorial function. The approximation is asymptotically correct for large n (written ~) which means that the ratio of n! and the approximation tends to 1 as n goes to infinity. Still, the approximation is quite accurate even for small values of n (we have two correct significant figures for n ≥ 7).

Stirling's formula is good for making estimates in combinatorics, and since combinatorial problems appear naturally in most branches of mathematics (particularly in probability theory) this makes Stirling's formula very useful.

Simple example:

The number of ways to partition a set with 3n members into three sets of size n is

(3n)!(n!)^-3 ~ 3^3n+1/2/2πn


In order to use the formula we should of course convince ourselves that it is correct.

As an informal justification for why there should be some formula like Stirling's we can begin by considering approximations to ln(n!). By comparing the sum ln(1) + ... + ln(n) with the integral ∫₀ⁿln(x)dx it is obvious that ln(n!) ~ n ln(n). This is the weak form of Stirling's formula. In order to derive the strong form we have to work a bit harder.

Proposition:

n! ~ (2πn)^1/2nⁿe^-n

Proof:

The proof consists of two parts (and a lot of dull equations). First we show that n! ~ Cn^n+1/2e^-n for some constant C, and then that the value of the constant is C = (2π)^1/2.

We know that for all positive x

1 - x + x² - x³ < 1/(1+x) < 1 - x + x²

Integrating this gives the inequality

x - x²/2 + x³/3 - x⁴/4 < ln(1+x) < x - x²/2 + x³/3

Now let h_n = ln(n!eⁿn^-n-1/2). Then

h_n - h_n+1 = ln(n!eⁿ(n+1)^n+3/2/(n+1)!eⁿ⁺¹n^n+1/2) = ln((n + 1/n)^n+1/2/e) = (n+1/2)ln(1 + 1/n) - 1

Using the estimation for ln(1+ 1/n) we obtain after some algebra

n^-2/12 - n^-3/12 - n^-4/8 < h_n - h_n+1 < n^-2/12 + n^-3/6

For n ≥ 2 we see that h_n - h_n+1 is positive, so the sequence h_n is a decreasing. Using the other inequality we can find a lower bound for the sequence by summing the series n^-2/12 + n^-3/6, which converges. By the completeness of the reals this implies that h_n tends to some finite limit A as n → ∞. Thus

n! ~ e^An^n+1/2e^-n

In order to complete the proof we need to pin down the value of A. An (admittedly somewhat unintuitive) way of doing this is to consider the sequence of integrals

I_n = ∫₀^π/2sinⁿθ dθ

I₀ = π/2, I₁ = 1 by direct integration, and integration by parts gives the recursion relation I_n = (n-1)I_n-2/n for n ≥ 2. Hence we obtain different expressions for I_n depending on whether n is even or odd:

I_2n = (2ⁿn!)^-2(2n)!π/2
I_2n+1 = (2ⁿn!)²/(2n+1)!

Using our proto-version n! ~ e^An^n+1/2e^-n of Stirling's formula we can write (after a lot of cancellations)

I_2n/I_2n+1 = (2n+1)((2n)!)²(2ⁿn!)^-4π/2 ~ (2n+1)(2/n)e^-2Aπ/2 → 2πe^-2A

But we also know that I_n is a decreasing sequence, so

1 < I_2n/I_2n+1 < I_2n-1/I_2n+1 = 1 +(2n)^-1 → 1

as n → ∞. Since limits are unique

2πe^-2A = 1
e^A = (2π)^1/2
n! ~ (2πn)^1/2nⁿe^-n

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

• Steepest Descent
• Asymptotic
• Exponentially Small