Modifications to Stirling's approximation

We consider a handful of useful approximations to the factorial or gamma function.

Stirling's approximation

We begin with a derivation of Stirling's approximation. It is easy to show, using integration by parts, that the following integral representation of the gamma function gives the factorial of integer values.

$n!={\int }_0^{\mathrm{\infty }}\mathrm{dt}{t}^{\mathrm{n}}{e}^{-t}=\Gamma (n+1)$

The Γ function is a continuous generalization of the integer-defined factorial function, with the definition such that $\Gamma \left(n\right)=(n-1)!$. This integral is also known as a Euler integral of the second kind. In order to derive Stirling's approximation, we first rewrite the integrand.

$n!={\int }_0^{\mathrm{\infty }}\mathrm{dt}{e}^{f\left(t\right)}={\int }_0^{\mathrm{\infty }}\mathrm{dt}{e}^{nlnt-t}$

For values of $t<1$, the exponential approaches zero because of the logarithm. For large values of $t$, the exponential also approaches zero. The integral is dominated by values of $t$ near the local maximum of $f\left(t\right)$. For $f\left(t\right)$ and its derivatives we find

$f\left(t\right)=nlnt-t$,    $f\text{'}\left(t\right)=n/t-1$   and    $f\text{'}\text{'}\left(t\right)=-n/t^2$.

The derivative tells us that the maximum is at $t=n$. Near the maximum we approximate the function in terms of $x=t-n$ as

$f\left(x\right)\approx nlnn-n-x^2/2n$,

which gives us an approximation to the integral as

$n!\approx n^{\mathrm{n}}{e}^{-n}{\int }_{-n}^{\mathrm{\infty }}\mathrm{dx}{e}^{-x^2/2n}.$

The usual application of this approximation is for large values of $n$, in which case we can treat the integral's lower limit as minus infinity without additional error. This gives the standard version of Stirling's approximation.

$n!\approx \sqrt{2\mathrm{<span\; class="pie">\pi </span>}n}{n}^{\mathrm{n}}{e}^{-n}$

Modified approximation

We begin with the same integral.

$n!={\int }_0^{\mathrm{\infty }}\mathrm{dt}{t}^{\mathrm{n}}{e}^{-t}$

This time, we rewrite the integral with a change of variables in terms of $x=lnt$. The result is

$n!={\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}\mathrm{dx}{e}^{(n+1)x-e^{\mathrm{x}}}$

We use the same treatment of this integral as the previous one. For $f\left(x\right)$ and its derivatives we find

$f\left(x\right)=(n+1)x-e^{\mathrm{x}}$,    $f\text{'}\left(x\right)=n+1-e^{\mathrm{x}}$   and    $f\text{'}\text{'}\left(x\right)=-e^{\mathrm{x}}$.

The derivative tells us that the maximum is at $e^{\mathrm{x}}=n+1$. Near the maximum we approximate the function in terms of $y=x-ln(n+1)$ as

$f\left(y\right)\approx (n+1)ln(n+1)-(n+1)-(n+1)y^2/2$ ,

with the corresponding approximation to the integral

$n!\approx {(n+1)}^{n+1}{e}^{-(n+1)}{\int }_{-\infty }^{\infty }\mathrm{dy}{e}^{-(n+1)y^2/2}.$

We arrive at our first modified version of Stirling's approximation:

$n!\approx \sqrt{2\mathrm{<span\; class="pie">\pi </span>}/(n+1)}{(n+1)}^{n+1}{e}^{-(n+1)}$

The following plot shows a comparison of these approximations with the gamma (factorial) function from 2 to 3. First, it is somewhat remarkable how well they both do at such small values of n, since derivations generally assume large values of n. The relevant factorials are obviously 2!=2 and 3!=6. Stirling's approximation gives 1.9190 and 5.8362, with errors of -4.0% and -2.7%. The second formula, expressed in terms of n+1, gives 1.9454 and 5.8765, for errors of -2.7% and -2.1%. There is a much more remarkable difference at n=0. The true value is 0!=1. Stirling's approximation gives 0. The modified formula gives $\sqrt{2\mathrm{<span\; class="pie">\pi </span>}}/e\approx \mathrm{0.92}$.

Comparison of gamma (factorial) function, Stirling's approximation, and the preceding equation from 2 to 3. The δ=0.5 label is a function also known as Burnside's formula.

Generalization between the two leading order approximations

If we compare Stirling's approximation with the modified formula, they can both be represented by the following equation

$n!\approx \sqrt{2\mathrm{<span\; class="pie">\pi </span>}}{(n+\delta )}^{n+1/2}{e}^{-(n+\delta )}$

Stirling's approximation is recovered with $\delta =0$. The modified formula follows for $\delta =1$. Interestingly, there is improved accuracy for intermediate values. The plot above includes the example of $\delta =\mathrm{0.5}$. Over this range, this is an improvement, with a value of 2.0333 at 2 and 6.0715 at 3, larger than the correct factorial values by 1.7% and 1.2%. Optimal values (in some regard) for δ can be found in Mortici. Below, we will show that these special values of δ cancel the 1/m² term in an asymptotic expansion.

Asymptotic expansion of Stirling's approximation

Now we consider the accuracy of Stirling's approximation in an asymptotic sense and derive a higher order term that makes it more accurate for large values of n. Consider the ratio of the true factorial to Stirling's approximation.

$r_{\mathrm{n}}=\frac{n!}{\sqrt{2\mathrm{<span\; class="pie">\pi </span>}n}{n}^{\mathrm{n}}{e}^{-\mathrm{n}}}$

If the function were precisely equal to the factorial, both this ratio and the ratio of sequential terms in this series would be unity. The ratio of sequential terms is

$\frac{r_{\mathrm{n}}}{r_{n-1}}=\frac{n!}{\sqrt{2\mathrm{<span\; class="pie">\pi </span>}n}{n}^{\mathrm{n}}{e}^{-n}}\frac{\sqrt{2\mathrm{<span\; class="pie">\pi </span>}(n-1)}{(n-1)}^{(n-1)}{e}^{-(n-1)}}{(n-1)!}$

This simplifies to

$\frac{r_{\mathrm{n}}}{r_{n-1}}=en\sqrt{\frac{(n-1)}{n}}\frac{{(n-1)}^{(n-1)}}{n^{\mathrm{n}}}=e\frac{{(n-1)}^{(n-1/2)}}{n^{(n-1/2)}}=e{(1-1/n)}^{n-1/2}$,

which we finally write as

$\frac{r_{\mathrm{n}}}{r_{n-1}}=e{(1-1/n)}^n{(1-1/n)}^{-1/2}=R$

If this ratio is one, then its logarithm will be zero.

$lnR=1+nln(1-1/n)-(1/2)ln(1-1/n)$

This expression is convenient because we can use the Taylor expansion of the logarithms to get the asymptotic dependence in powers of 1/n.

$nln(1-1/n)=n[-1/n-1/(2n^2)-1/(3n^3)-\mathrm{<span\; class="pie">\dots </span>}]$,

with the ellipsis denoting neglect of higher powers of 1/n. Similarly, the second logarithm can be written

$(1/2)ln(1-1/n)=(1/2)[-1/n-1/(2n^2)-1/(3n^3)-\mathrm{<span\; class="pie">\dots </span>}]$.

If we combine the terms, then to second order in 1/n, we find

$lnR=1-n/n-n/\left(2n^2\right)-n/\left(3n^3\right)+(1/2)1/n+(1/2)1/\left(2n^2\right)+\mathrm{<span\; class="pie">\dots </span>}$

Simplifying the combined terms, we find

$lnR=-1/\left(\mathrm{12}{n}^2\right)+$ (terms proportional to $1/n^3$)

It is not difficult to show, by following the same procedure, that the following modification to Stirling's approximation would cancel this term in $lnR$:

$n!\approx \sqrt{2\mathrm{<span\; class="pie">\pi </span>}n}{n}^{\mathrm{n}}{e}^{-n+1/\left(\mathrm{12}n\right)}$

Asymptotic expansion of modified formula

The same asymptotic expansion can be applied to the modified formula. When we do this, we find

$n!\approx \sqrt{2\mathrm{<span\; class="pie">\pi </span>}/(n+1)}{(n+1)}^{n+1}{e}^{-(n+1)+1/\left(\mathrm{12}\right(n+1\left)\right)}$

This form makes the gamma function notation too convenient, so finally we write

$\Gamma \left(x\right)=(x-1)!\approx \sqrt{2\mathrm{<span\; class="pie">\pi </span>}}{x}^{\mathrm{x}-1/2}{e}^{-x+1/\left(\mathrm{12}x\right)}$

The gamma function is often written with an argument of z to denote that it is defined throughout the complex plane. Here it is written as a function of x, because the approximation is valid for the factorial of real numbers greater than zero (corresponding to gamma function arguments x ≥ 1).

Generalized asymptotic expansion

Let's return to the generalized formula

$n!\approx \sqrt{2\mathrm{<span\; class="pie">\pi </span>}}{(n+\delta )}^{n+1/2}{e}^{-(n+\delta )}$

We can perform another asymptotic study on this functional form. Along the way, we define m=n+δ. Then for the same power series expansion performed on the standard Stirling's approximation, we find the first two non-zero terms to be

$lnR=1/m^2(-1/\mathrm{12}+\delta /2-{\delta }^2/2)+1/m^3(-1/\mathrm{12}+\delta /3-{\delta }^3/3)+$ (terms proportional to $1/m^4$)

The first term is zero for values $\delta =(1/2)(1\pm 1/\sqrt{3})$, which happen to be the values given in Mortici. So these optimal values, as given, cancel the second order term. However, these values result in a non-zero 1/m³ term. Interestingly there are three values of δ which result in these two terms being equal: 0, 1/2 and 1. For these values, the 1/m² term as given above also cancels the 1/m³ term.

Euler-Maclaurin formula result

We can also apply the Euler-Maclaurin formula directly to the logarithm of n!. We set the integration constant to half the logarithm of 2π based on previous considerations. Including terms through the ninth derivative of the logarithm, this yields

$n!\approx \sqrt{2\mathrm{<span\; class="pie">\pi </span>}n}{n}^{\mathrm{n}}{e}^{-n+f(n)}$

$f(n)=\frac{1}{\mathrm{12}n}-\frac{1}{\mathrm{360}{n}^3}+\frac{1}{\mathrm{1260}{n}^5}-\frac{1}{\mathrm{1680}{n}^7}+\frac{1}{\mathrm{1188}{n}^9}$

This is identical to 6.1.41 in Handbook of Mathematical Functions, with one additional term included.

Comparison and summary

The advantage in the modified formula is due to the form the integral takes before application of the approximation. The limits extend to both plus and minus infinity after subsitution, plus the functional form of the exponential is better suited to the approximation. It actually even works for values of x < 1 in the argument of the gamma function, but it fails as the gamma function approaches x=0.

Behavior of the factorial function and approximations from 0 to 2. The correct value of 0! is 1. Stirling's approximation gives 0. The higher order (H.O.) version of Stirling's approximation diverges near zero, but agrees much better than the original version as x approaches 0.5. The modified formula is a decent approximation with an offset less than 8% throughout this range. The H.O. version of the modified formula is nearly indistinguishable from the correct value throughout this range.

We now show a precise comparison between some of the approximations given here and some very good ones. The one from Luschny looks like the Burnside formula, but with four correction terms in continued fraction form. The one labelled δ=0.7887 is one of the versions that cancel the 1/m² expansion term. The relatively simple modified HO form derived here is practically identical to 6.1.37 in Handbook of Mathematical Functions including terms through 1/(288z²) in square brackets. It is the best at the smallest values, then the Luschny form is best until about x=4. The Euler-Maclaurin formula result ('Euler') is the best above x≈4 and rapidly approaches floating point precision above x=10. It is amusing that the simple form approximates the factorial function well to x values as low as -0.5 (gamma function argument of 0.5).

Relative errors in five different approximations to the factorial function.

[1] "An ultimate extremely accurate formula for approximation of the factorial function," by Cristinel Mortici, in Archiv der Mathematik 93 (2009).
[2] "Approximation Formulas for the Factorial Function n!," by Peter Luschny
[3] "Handbook of Mathematical Functions," Abramowitz and Stegun, Dover Publications

• Mathematics and physics
• Ski Sickness
• Patient history - chronological index