Stirling's formula: proof?

Question

Suppose we want to show that $$ n! \sim \sqrt{2 \pi} n^{n+(1/2)}e^{-n}$$

Instead we could show that $$\lim_{n \to \infty} \frac{n!}{n^{n+(1/2)}e^{-n}} = C$$ where $C$ is a constant. Maybe $C = \sqrt{2 \pi}$.

What is a good way of doing this? Could we use L'Hopital's Rule? Or maybe take the log of both sides (e.g., compute the limit of the log of the quantity)? So for example do the following $$\lim_{n \to \infty} \log \left[\frac{n!}{n^{n+(1/2)}e^{-n}} \right] = \log C$$

Answer 1

A proof I found a while ago entirely relies on creative telescoping. Since $\frac{1}{n^2}-\frac{1}{n(n+1)}=\frac{1}{n^2(n+1)}$,

$$\begin{eqnarray*} \sum_{n\geq m}\frac{1}{n^2}&=&\sum_{n\geq m}\left(\frac{1}{n}-\frac{1}{(n+1)}\right)+\frac{1}{2}\sum_{n\geq m}\left(\frac{1}{n^2}-\frac{1}{(n+1)^2}\right)\\&+&\frac{1}{6}\sum_{n\geq m}\left(\frac{1}{n^3}-\frac{1}{(n+1)^3}\right)-\frac{1}{6}\sum_{n\geq m}\frac{1}{n^3(n+1)^3}\tag{1}\end{eqnarray*} $$ hence, by the series representation for $\psi(z)=\frac{d}{dz}\log\Gamma(z)$ (where $\Gamma(z)$ is the analytic continuation of $\int_{0}^{+\infty}t^{z-1}e^{-t}\,dt$, defined for $\text{Re}(z)>0$): $$ \psi'(m)=\sum_{n\geq m}\frac{1}{n^2}\leq \frac{1}{m}+\frac{1}{2m^2}+\frac{1}{6m^3}\tag{2}$$ and in a similar fashion: $$ \psi'(m) \geq \frac{1}{m}+\frac{1}{2m^2}+\frac{1}{6m^3}-\frac{1}{30m^5}.\tag{3}$$ Integrating twice, we have that $\log\Gamma(m)$ behaves like: $$ \log\Gamma(m)\approx\left(m-\frac{1}{2}\right)\log(m)-\color{red}{\alpha} m+\color{blue}{\beta}+\frac{1}{12m}\tag{4}$$ where $\color{red}{\alpha=1}$ follows from $\log\Gamma(m+1)-\log\Gamma(m)=\log m$.

That gives Stirling's inequality up to a multiplicative constant.

$\color{blue}{\beta=\log\sqrt{2\pi}}$ then follows from Legendre's duplication formula and the well-known identity:

$$ \Gamma\left(\frac{1}{2}\right)=2 \int_{0}^{+\infty}e^{-x^2}\,dx = \sqrt{\pi}.\tag{5}$$

---

Addendum: if we apply creative telescoping like in the second part of this answer, i.e. by noticing that $k(x)=\frac{60x^2-60x+31}{60x^3-90x^2+66x-18}$ gives $k(x)-k(x+1)=\frac{1}{x^2}+O\left(\frac{1}{x^8}\right)$, we arrive at

$$\begin{eqnarray*} m!&\approx& 2^{\frac{37-32m}{42}}e^{\frac{1}{84} \left(42-\sqrt{35} \pi -84 m+2 \sqrt{35} \arctan\left[\sqrt{\frac{5}{7}} (2m-1)\right]\right)} \\ &\cdot&\sqrt{\pi}\, m\,(2m-1)^{\frac{8}{21}(2m-1)}\left(m^2-m+\frac{3}{5}\right)^{\frac{5}{84}(2m-1)}\tag{6}\end{eqnarray*} $$ that is much more accurate than the "usual" Stirling's inequality, but also way less "practical".
However, it might be fun to plug in different values of $m$ in $(6)$ to derive bizarre approximate identities involving $e,\pi,\sqrt{\pi}$ and values of the arctangent function, like

$$ \sqrt{\pi}\,\exp\left[-\frac{1}{42}\left(147+\sqrt{35} \arctan\frac{1}{\sqrt{35}}\right)\right]\approx 2^{\frac{19}{6}}3^{\frac{1}{6}}5^{\frac{5}{12}} 7^{-\frac{37}{12}}.\tag{7}$$

Answer 2

Inspired by the two references below, here's a short proof stripped of motivation and details.

For $t>0$, define

$$ g_t(y) = \begin{cases} \displaystyle \left(1+\frac{y}{\sqrt{t}}\right)^{\!t} \,e^{-y\sqrt{t}} & \text{if } y>-\sqrt{t}, \\ 0 & \text{otherwise}. \end{cases} $$

It is not hard to show that $t\mapsto g_t(y)$ is decreasing for $y\geq 0$ and increasing for $y\leq 0$. The limit is $g_\infty(y)=\exp(-y^2/2)$, so dominated convergence gives

$$ \lim_{t\to\infty}\int_{-\infty}^\infty g_t(y)\,dy = \int_{-\infty}^\infty \exp(-y^2/2)\,dy = \sqrt{2\pi}. $$

Use the change of variables $x=y\sqrt{t}+t$ to get

$$ t! = \int_0^\infty x^t e^{-x}\,dx = \left(\frac{t}{e}\right)^t \sqrt{t} \int_{-\infty}^\infty g_t(y)\,dy. $$

References:

[1] J.M. Patin, A very short proof of Stirling's formula, American Mathematical Monthly 96 (1989), 41-42.

[2] Reinhard Michel, The $(n+1)$th proof of Stirling's formula, American Mathematical Monthly 115 (2008), 844-845.

Answer 3

Here is a derivation of Stirling's Formula as given by

$$\bbox[5px,border:2px solid #C0A000]{\sqrt{2\pi n}\left(\frac ne\right)^n\left(1+\frac{1}{12(n+1)}\right) \le n!\le \sqrt{2\pi n}\left(\frac ne\right)^n\left(1+\frac{1}{12(n-2)}\right)}\tag 1$$

In the ensuing development, we will use $(i)$ the formula for the trapezoidal rule for twice differentiable functions, $(ii)$ Wallis's Formula for $\pi$, and $(iii)$ standard inequalities.

---

USING THE TRAPEZOIDAL RULE:

To prove $(1)$, we begin with the trapezoidal rule formula

$$\int_a^b f(x)\,dx=\frac12\left(f(a)+f(b)\right)(b-a)-\frac1{12}f''(\xi)(b-a)^3\tag 2$$

for which $f''(x)$ exists for $x\in[a,b]$ and $\xi\in (a,b)$. Letting $f(x)=\log(x)$ in $(2)$ shows that

$$\begin{align} \int_1^n \log(x)\,dx&=\sum_{m=1}^{n-1} \int_m^{m+1}\log(x)\,dx\\\\ &=\sum_{m=1}^{n-1} \left(\frac12 (\log(m)+\log(m+1))+\frac1{12}\frac{1}{\xi_m^2}\right)\tag 3 \end{align}$$

where $m<\xi_m<m+1$. Denoting $\frac1{12}\sum_{m=1}^{n-1}\frac{1}{\xi_m^2}$ by $C_n$, we see that $\lim_{n\to \infty}C_n=C<\infty$.

Now, we can rewrite $(3)$ as

$$\log(n!)=\log\left(e^{1-C_n}\sqrt{n}\left(\frac{n}{e}\right)^n\right)$$

which becomes

$$\bbox[5px,border:2px solid #C0A000]{n!=e^{1-C_n}\sqrt{n}\left(\frac{n}{e}\right)^n }\tag 4$$

---

USING WALLIS'S PROUCT FORMULA:

To evaluate the limit $\lim_{n\to \infty}e^{1-C_n}$, we rely on Wallis's Product for $\pi$ as given by

$$\frac{\pi}{2}=\lim_{n\to \infty}\frac{2^{4n}(n!)^4}{(2n+1)((2n)!)^2}\tag 5$$

Using $(4)$ in $(5)$ shows that

$$\frac{\pi}{2}=\lim_{n\to \infty}\frac{n^2}{(2n)(2n+1)}\left(e^{1-2C_n+C_{2n}}\right)^2$$

from which we find that

$$\bbox[5px,border:2px solid #C0A000]{\lim_{n\to \infty}e^{1-C_n}=\sqrt{2\pi}}\tag 6$$

---

PUTTING IT TOGETHER TO ARRIVE AT STIRLING'S FORMULA

We are now able to write the factorial

$$n!=\sqrt{2\pi n}\left(\frac ne\right)^n e^{C_n-C} \tag 7$$

Since $C_n\to C$, then we have Stirling's Formula

$$\bbox[5px,border:2px solid #C0A000]{n!\sim \sqrt{2\pi n}\left(\frac ne\right)^n }$$

as was to be shown!

---

DERIVING BOUNDS FOR $\displaystyle n!$:

It is easy to see that since $\frac{1}{(m+1)^2}<\frac1{\xi_m^2}<\frac{1}{m^2}$, then

$$C-C_n<\frac1{12}\int_{n-1}^\infty \frac{1}{x^2}\,dx=\frac{1}{12(n-1)}$$

and

$$C-C_n>\frac1{12}\int_{n+1}^\infty \frac{1}{x^2}\,dx=\frac{1}{12(n+1)}$$

Next, for $x<1/2$ the Mean-Value Theorem guarantees that

$$1+x<e^x<1+x+x^2$$

Hence, we see that for $n\ge 3$

$$1+\frac{1}{12(n+1)}<e^{C-C_n}<1+\frac{1}{12(n-1)}+\frac{1}{144(n-1)^2}<1+\frac{1}{12(n-2)}\tag 8$$

Finally, using $(8)$ in $(7)$ yields the coveted bounds in $(1)$

$$\bbox[5px,border:2px solid #C0A000]{\sqrt{2\pi n}\left(\frac ne\right)^n\left(1+\frac{1}{12(n+1)}\right) \le n!\le \sqrt{2\pi n}\left(\frac ne\right)^n\left(1+\frac{1}{12(n-2)}\right)}$$

Answer 4

The way I've usually shown Stirling's Asymptotic Approximation starts with the formula $$ \begin{align} n! &=\int_0^\infty x^ne^{-x}\;\mathrm{d}x\\ &=\int_0^\infty e^{-x+n\log(x)}\;\mathrm{d}x\\ &=\int_0^\infty e^{-nx+n\log(nx)}n\;\mathrm{d}x\\ &=n^{n+1}\int_0^\infty e^{-nx+n\log(x)}\;\mathrm{d}x\\ &=n^{n+1}e^{-n}\int_{-1}^\infty e^{-nx+n\log(1+x)}\;\mathrm{d}x\\ &=n^{n+1}e^{-n}\int_{-\infty}^\infty e^{-nu^2/2}x'\;\mathrm{d}u\tag{1} \end{align} $$ Where $u^2/2=x-\log(1+x)$. That is, $u(x)=x\sqrt{\frac{x-\log(1+x)}{x^2/2}}$, which is analytic near $x=0$: the singularity of $\frac{x-\log(1+x)}{x^2/2}$ at $x=0$ is removable, and since $\frac{x-\log(1+x)}{x^2/2}$ is near $1$ when $x$ is near $0$, $\sqrt{\frac{x-\log(1+x)}{x^2/2}}$ is analytic near $x=0$. Since $u(0)=0$ and $u'(0)=1$, the Lagrange Inversion Theorem says that $x$ is an analytic function of $u$ near $u=0$, with $x(0)=0$ and $x'(0)=1$.

Note that since $u^2/2=x-\log(1+x)$, we have $$ u(1+x)=xx'\tag{2} $$ from which it follows that $\lim\limits_{u\to+\infty}\dfrac{x'(u)}{u}=1$ and $\lim\limits_{u\to-\infty}\dfrac{x'(u)}{ue^{-u^2/2}}=-1$. That is, $x'(u)=O(u)$ for $|u|$ near $\infty$.

Because $x$ is an analytic function with $x(0)=0$ and $x'(0)=1$, $$ x=\sum_{k=1}^\infty a_ku^k\tag{3} $$ where $a_1=1$. Then, looking at the coefficient of $u^n$ in $(2)$, we get $$ \begin{align} a_{n-1} &=\sum_{k=1}^na_{n-k+1}ka_k\\ &=\sum_{k=1}^na_k(n-k+1)a_{n-k+1}\\ &=\frac{n+1}{2}\sum_{k=1}^na_ka_{n-k+1}\tag{4} \end{align} $$ Therefore, we get the recursion $$ a_n=\frac{a_{n-1}}{n+1}-\frac{1}{2}\sum_{k=2}^{n-1}a_ka_{n-k+1}\tag{5} $$ Thus, we get the beginning of the power series for $x(u)$ to be $$ x=u+\tfrac{1}{3}u^2+\tfrac{1}{36}u^3-\tfrac{1}{270}u^4+\tfrac{1}{4320}u^5+\tfrac{1}{17010}u^6-\tfrac{139}{5443200}u^7+\tfrac{1}{204120}u^8+\dots\tag{6} $$ Integration by parts yields the following two identities: $$ \int_{-\infty}^\infty e^{-nu^2}|u|^{2k+1}\;\mathrm{d}u=\frac{k!}{n^{k+1}}\tag{7a} $$ and $$ \int_{-\infty}^\infty e^{-nu^2}u^{2k}\;\mathrm{d}u=\frac{(2k-1)!!}{2^kn^{k+1/2}}\sqrt{\pi}\tag{7b} $$ Furthermore, we have the following tail estimates: $$ \begin{align} \int_{|u|>a}e^{-nu^2}\;\mathrm{d}u &=\int_a^\infty e^{-nu^2}u^{-1}\;\mathrm{d}u^2\\ &=\int_{a^2}^\infty e^{-nu}u^{-1/2}\;\mathrm{d}u\\ &\le\frac{1}{a}\int_{a^2}^\infty e^{-nu}\;\mathrm{d}u\\ &=\frac{1}{na}\int_{na^2}^\infty e^{-u}\;\mathrm{d}u\\ &=\frac{1}{na}e^{-na^2}\tag{8a} \end{align} $$ and for $k\ge1$, $$ \begin{align} \int_{|u|>a}e^{-nu^2}|u|^{k}\;\mathrm{d}u &=\int_a^\infty e^{-nu^2}u^{k-1}\;\mathrm{d}u^2\\ &=\int_{a^2}^\infty e^{-nu}u^{(k-1)/2}\;\mathrm{d}u\\ &=n^{-(k+1)/2}\int_{na^2}^\infty e^{-u}u^{(k-1)/2}\;\mathrm{d}u\\ &=n^{-(k+1)/2}\int_{na^2}^\infty e^{-u/2}e^{-u/2}u^{(k-1)/2}\;\mathrm{d}u\\ &\le n^{-(k+1)/2}\int_{na^2}^\infty e^{-u/2}\left(\frac{k-1}{e}\right)^{(k-1)/2}\;\mathrm{d}u\\ &=\frac{2}{n}\left(\frac{k-1}{ne}\right)^{(k-1)/2}e^{-na^2/2}\tag{8b} \end{align} $$ The tail estimates in $(8)$ are $O\left(\frac{1}{n}e^{-na^2/2}\right)$ for all $k\ge0$. That is, they decay faster than any power of $n$.

Define $$ f_k(u)=x'(u)-a_1-2a_2u-\dots-2ka_{2k}u^{2k-1}\tag{9} $$ Since $x(u)$ is analytic near $u=0$, $f_k(u)=O\left(u^{2k}\right)$ near $u=0$. Since $x'(u)=O(u)$ for $|u|$ near $\infty$, $f_k(u)=O\left(u^{2k-1}\right)$ for $|u|$ near $\infty$. Therefore, $$ \begin{align} \int_{-\infty}^\infty e^{-nu^2/2}f_k(u)\;\mathrm{d}u &=\int_{|u|< a}e^{-nu^2/2}f_k(u)\;\mathrm{d}u + \int_{|u|> a}e^{-nu^2/2}f_k(u)\;\mathrm{d}u\\ &\le\int_{|u|< a}e^{-nu^2/2}C_1u^{2k}\;\mathrm{d}u + \int_{|u|> a}e^{-nu^2/2}C_2|u|^{2k-1}\;\mathrm{d}u\\ &=O\left(\frac{1}{n^{k+1/2}}\right)\tag{10} \end{align} $$ Therefore, using $\mathrm{(7b)}$ and $(10)$, we get $$ \begin{align} n! &=n^{n+1}e^{-n}\int_{-\infty}^\infty e^{-nu^2/2}\left(1+\tfrac{1}{12}u^2+\tfrac{1}{864}u^4-\tfrac{139}{777600}u^6+f_4(u)\right)\;\mathrm{d}u\\ &+\;n^{n+1}e^{-n}\int_{-\infty}^\infty e^{-nu^2/2}\left(\tfrac{2}{3}u-\tfrac{2}{135}u^3+\tfrac{1}{2835}u^5+\tfrac{1}{25515}u^7\right)\;\mathrm{d}u\\ &=\sqrt{2n}\;n^ne^{-n}\left(\int_{-\infty}^\infty e^{-u^2}\left(1+\tfrac{1}{6n}u^2+\tfrac{1}{216n^2}u^4-\tfrac{139}{97200n^3}u^6\right)\;\mathrm{d}u+O\left(\frac{1}{n^4}\right)\right)\\ &=\sqrt{2\pi n}\;n^ne^{-n}\left(1+\frac{1}{12n}+\frac{1}{288n^2}-\frac{139}{51840n^3}+O\left(\frac{1}{n^4}\right)\right)\tag{11} \end{align} $$

Answer 5

Instead of a proof, how about a string of hints? This comes from Maxwell Rosenlicht's Introduction to Analysis (a great little easy-to-read text which is dirt cheap -- it's a Dover paperback).

• Chapter VI: Problem #22 Show that for $n=1,2,3,\dots$ we have that $$ 1+\frac{1}{2}+\frac{1}{3}+\cdots+\frac{1}{n}-\ln(n)$$ is positive and it decreases as $n$ increases. Thus this converges to a number between 0 and 1 (Euler's constant).

• Chapter VII: Problem #39 For $n=0,1,2,\dots$ let $I_n=\int_0^{\pi/2} \sin^n(x)\,dx$. Show that

(a) $\frac{d}{dx}\left(\cos(x)\sin^{n-1}(x)\right) = (n-1)\sin^{n-2}(x)-n\sin^n(x)$

(b) $I_n = \frac{n-1}{n}I_{n-2}$ if $n \geq 2$

(c) $I_{2n} = \frac{1\cdot 3\cdot 5\cdots (2n-1)}{2\cdot 4\cdot 6\cdots (2n)}\frac{\pi}{2}$ and $I_{2n+1} = \frac{2\cdot 4\cdot 6\cdots (2n)}{3\cdot 5\cdot 7\cdots (2n+1)}$ for $n=1,2,3,\dots$

(d) $I_0, I_1, I_2, \dots$ is a decreasing sequence having the limit zero and $$ \lim_{n \to \infty} \frac{I_{2n+1}}{I_{2n}}=1$$

(e) Wallis' product: $$\lim_{n \to \infty} \frac{2\cdot 2\cdot 4\cdot 4\cdot 6\cdot 6 \cdots (2n)\cdot (2n)}{1\cdot 3\cdot 3\cdot 5\cdot 5\cdot 7 \cdots (2n-1)\cdot (2n+1)}=\frac{\pi}{2}$$

• Chapter VII: Problem #40

(a) Show that if $f:\{ x\in \mathbb{R}\;|\; x\geq 1\} \to \mathbb{R}$ is continuous, then $$ \sum\limits_{i=1}^n f(i) = \int_1^{n+1} f(x)\,dx + \sum\limits_{i=1}^n\left(f(i)-\int_i^{i+1} f(x)\,dx\right)$$

(b) Show that if $i>1$, then $\ln(i)-\int_i^{i+1}\,dx$ differs from $-1/2i$ by less than $1/6i^2$. [Hint: Work out the integral using the Taylor series for $\ln(1+x)$ at the point $0$.]

(c) Use part (a) with $f=\ln$, part (b), and Problem #22 from Chapter VI to prove that the following limit exists: $$\lim_{n \to \infty} \left(\ln(n!)-\left(n+\frac{1}{2}\right)\ln(n)+n\right)$$

(d) Use part (e) of Problem #39 to compute the above limit, then obtaining: $$ \lim_{n\to \infty} \frac{n!}{n^ne^{-n}\sqrt{2\pi n}}=1$$ (i.e. Stirling's Formula)

Answer 6

Here are a couple of "proofs" of Stirling's formula. They are quite elegant (in my opinion), but not rigorous. On could write down a real proof from these, but as they rely on some hidden machinery, the result would be quite heavy.

1) A probabilistic non-proof

We start from the expression $e^{-n} n^n/n!$, of which we want to find an equivalent. Let us fix $n$, and let $Y$ be a random variable with a Poisson distribution of parameter $n$. By definition, for any integer $k$, we have $\mathbb{P} (Y=k) = e^{-n} n^k/k!$. If we take $k=n$, we get $\mathbb{P} (Y=n) = e^{-n} n^n/n!$. The sum of $n$ independent random variables with a Poisson distribution of parameter $1$ has a Poisson distribution of parameter $n$; so let us take a sequence $(X_k)$ of i.i.d. random variables with a Poisson distribution of parameter $1$. Note that $\mathbb{E} (X_0) = 1$. We have:

$$\mathbb{P} \left( \sum_{k=0}^{n-1} (X_k - \mathbb{E} (X_k)) = 0 \right) = \frac{e^{-n} n^n}{n!}.$$

In other words, $e^{-n} n^n/n!$ is the probability that a centered random walk with Poissonnian steps of parameter $1$ is in $0$ at time $n$. We have tools to estimates such quantities, namely local central limit theorems. They assert that:

$$\frac{e^{-n} n^n}{n!} = \mathbb{P} \left( \sum_{k=0}^{n-1} (X_k - \mathbb{E} (X_k)) = 0 \right) \sim \frac{1}{\sqrt{2 \pi n \text{Var} (X_0)}},$$

a formula which is closely liked with Gauss integral and diffusion processes. Since the variance of $X_0$ is $1$, we get:

$$n! \sim \sqrt{2 \pi n} n^n e^{-n}.$$

The catch is of course that the local central limit theorems are in no way elementary results (except for the simple random walks, and that is if you already know Stirling's formula...). The methods I know to prove such results involve Tauberian theorems and residue analysis. In some way, this probabilistic stuff is a way to disguise more classical approaches (in my defense, if all you have is a hammer, everything looks like a nail).

I think one could get higher order terms for Stirling's formula by computing more precise asymptotics for Green's function in $0$, which requires the knowledge of higher moments for the Poisson distribution. Note that the generating function for a Poisson distribution of parameter $1$ is:

$$\mathbb{E} (e^{t X_0}) = e^{e^t-1},$$

and this "exponential of exponential" will appear again in a moment.

2) A generating functions non-proof

If you want to apply analytic methods to problems related to sequences, generating functions are a very useful tool. Alas, the series $\sum_{n \geq 0} n! z^n$ is not convergent for non-zero values of $z$. Instead, we shall work with:

$$e^z = \sum_{n \geq 0} \frac{z^n}{n!};$$

we are lucky, as this generating function is well-known. Let $\gamma$ be a simple loop around $0$ in the complex plane, oriented counter-clockwise. Let us fix some non-negative integer $n$. By Cauchy's formula,

$$\frac{1}{n!} = \frac{1}{n!} \frac{\text{d} e^z}{\text{d} z}_{|_0} = \frac{1}{2 \pi i} \oint_\gamma \frac{e^z}{z^{n+1}} \text{d} z.$$

We choose for $\gamma$ the circle of radius $n$ around $0$, with its natural parametrization $z (t) = n e^{it}$:

$$\frac{1}{n!} = \frac{1}{2 \pi i} \int_{- \pi}^{\pi} \frac{e^{n e^{it}}}{n^{n+1} e^{(n+1)it}} nie^{it} \text{d} t = \frac{e^n}{2 \pi n^n} \oint_{- \pi}^{\pi} e^{n (e^{it}-it-1)} \text{d} t = \frac{e^n}{2 \pi \sqrt{n} n^n} \int_{- \pi \sqrt{n}}^{\pi \sqrt{n}} e^{n \left(e^{\frac{i\theta}{\sqrt{n}}}-\frac{i\theta}{\sqrt{n}}-1\right)} \text{d} \theta,$$

where $\theta =t \sqrt{n}$. Hitherto, we have an exact formula; note that we meet again the "exponential of exponential". Now comes the leap of faith. For $x$ close to $0$, the value of $e^x-x-1$ is roughly $x^2/2$. Moreover, the bounds of the integral get close to $- \infty$ and $ \infty$. Hence, for large $n$, we have:

$$\frac{1}{n!} \sim \frac{e^n}{2 \pi \sqrt{n} n^n} \int_{- \infty}^{+ \infty} e^{\frac{n}{2} \left(\frac{i\theta}{\sqrt{n}}\right)^2} \text{d} \theta = \frac{e^n}{2 \pi \sqrt{n} n^n} \int_{- \infty}^{+ \infty} e^{-\frac{\theta^2}{2}} \text{d} \theta = \frac{e^n}{\sqrt{2 \pi n} n^n}.$$

Of course, it is not at all trivial to prove that the equivalents we took are rigorous. Indeed, if one apply this method to bad generating functions (e.g. $(1-z)^{-1}$), they can get false results. However, this can be done for some admissible functions, and the exponential is one of them.

I have learnt this method thanks to Don Zagier. It is also explained in Generatingfunctionology, Chapter $5$ (III), where the author credits Hayman. The original reference seems to be A generalisation of Stirling's formula (Hayman, 1956), but I can't read it now.

One of the advantage of this method is that it becomes very easy to get the next terms in the asymptotic expansion of $n!$. You just have to develop further the function $e^x-x-1$ at $0$. Another advantage is that it is quite general, as it can be applied to many other sequences.

Answer 7

Depending on one's preferences, one might care that it is possible to understand Stirling's formula (perhaps really due to Laplace?) in a usefully more general context, namely, as an example of a "Laplace's method" or "stationary phase" sort of treatment of asymptotics of integrals.

This is available on-line on my page http://www.math.umn.edu/~garrett/m/v/, with the file being http://www.math.umn.edu/~garrett/m/v/basic_asymptotics.pdf

One might object to certain very-classical treatments which make Gamma appear as a singular thing. While I agree that it is of singular importance in applications within and without mathematics, the means of establishing its asymptotics are not.

Answer 8

The best derivation of Stirling's approximation I have seen starts from Euler's Summation Formula, valid for integers $a\leq b$ and $m\geq 1$ $$ \sum_{a\leq k<b}f(k) = \int_a^b f(x)\,dx +\sum_{k=1}^m \frac{B_k}{k!} f^{(k-1)}(b) -\sum_{k=1}^m \frac{B_k}{k!} f^{(k-1)}(a) \\+ (-1)^{m+1} \int_a^b \frac{B_m}{m!}f^{(m)}(x)\,dx$$

Here, $f^{(n)}(x)$ refers to the $n$-th derivative of $f$ evaluated at $x$, and $B_k$ is the $k$-th Bernoulli number.

The "remainder term" $(-1)^{m+1} \int_a^b \frac{B_m}{m!}f^{(m)}(x)\,dx$ lies between $0$ and the first discarded term.

Applying Euler's Summation Formula to $f(x) = \ln(x)$, and using $a=1, b=n$, we obtain $$ \sum_{1\leq k<n} \ln k = n\ln n - n +\sigma - \frac{\ln n}{2} + \sum_{k=1}^m \frac{B_{2k}}{2k(2k-1)n^{2k-1}} + \varphi_{m,n}\frac{B_{2m+2}}{(2m+2)(2m+1)n^{2m+1}} $$ where $\varphi_{m,n}$ depends on the $n$ and the number of terms $m$, but is always between zero and one. $\sigma$ is a constant which is not determined by application of the formula. The best method I have seen for determining that $\sigma = \ln(\sqrt{2\pi})$ is in Brylan Schmuland's answer above, although in Concrete Mathematics the derive it by applying the summation formula to $\sum\binom{2n}{k}$.

Adding $\ln n$ to both sides, we get $$ \ln n! = n\ln n - n + \frac{\ln n}{2} + \ln \sqrt{2\pi} + \frac1{12n}-\frac1{360n^3} + \frac\epsilon{n^5} $$ with an error term having $0 \leq \epsilon \leq \frac1{1260}$.

Lastly, exponentiate and series expand the terms of order $\frac1n$ or less to get $$ n! = \sqrt{2\pi n}\left(\frac{n}{e}\right)^n \left( 1 + \frac1{12n} + \frac1{288n^2} - \frac{139}{51840n^3} -\frac{571}2488320 n^4 + \frac{\rho}{n^5}\right) $$ with $-0.00001 < \rho < 0.00078$

At this level of approximation, the approximation, rounded to the nearest integer, is exact for up to $n=11$ (of course, the relative error in the approximation falls as $\frac1{n^5}$ with a small coefficient, and so is extremely small.

Answer 9

If you're familiar with

$$\mathop {\lim }\limits_{n \to \infty } \frac{{\left( {2n} \right)!!}}{{\left( {2n - 1} \right)!!}}\frac{1}{{\sqrt n }} = \sqrt \pi $$

Then you can use

$$\eqalign{ & \mathop {\lim }\limits_{n \to \infty } \frac{{\left( {2n} \right)!{!^2}}}{{\left( {2n} \right)!}}\frac{1}{{\sqrt n }} = \sqrt \pi \cr & \mathop {\lim }\limits_{n \to \infty } \frac{{{2^{2n}}{{\left( {n!} \right)}^2}}}{{\left( {2n} \right)!}}\frac{1}{{\sqrt n }} = \sqrt \pi \cr} $$

Now you can check that

$$\alpha = \mathop {\lim }\limits_{n \to \infty } \frac{{n!{e^n}}}{{{n^n}\sqrt n }} = \mathop {\lim }\limits_{n \to \infty } \frac{{\left( {2n} \right)!{e^{2n}}}}{{{{\left( {2n} \right)}^{2n}}\sqrt {2n} }}$$

exists. Then square the first expression and divide by the latter to get

$$\alpha = \mathop {\lim }\limits_{n \to \infty } \frac{{{{\left( {n!} \right)}^2}{e^{2n}}}}{{{n^{2n}}n}}\frac{{{{\left( {2n} \right)}^{2n}}\sqrt {2n} }}{{\left( {2n} \right)!{e^{2n}}}} = \mathop {\lim }\limits_{n \to \infty } \frac{{{{\left( {n!} \right)}^2}{2^{2n}}\sqrt 2 }}{{\left( {2n} \right)!\sqrt n }} = \sqrt {2\pi } $$

Thus you have that

$$\mathop {\lim }\limits_{n \to \infty } \frac{{n!{e^n}}}{{{n^n}\sqrt {2n} }} = \sqrt \pi $$

or

$$n! \sim {n^n}{e^{ - n}}\sqrt {2\pi n} $$

Answer 10

This answer derives Stirling's Asymptotic Approximation by writing the log of a factorial as a Riemann-Stieltjes Integral.

---

Theorem $\bf{1}$: (Asymptotic Approximation of the Central Binomial Coefficients) $$ \lim_{n\to\infty}\frac{\sqrt{\pi n}}{4^n}\binom{2n}{n}=1\tag1 $$

Note that Theorem $1$ is equivalent to Wallis' Product for $\pi$.

Proof: The orthogonality of $e^{ikx}$ on $[0,2\pi]$ gives $$ \begin{align} \int_0^{2\pi}\cos^{2n}(x)\,\mathrm{d}x &=\frac1{4^n}\sum_{k=0}^{2n}\binom{2n}{k}\int_0^{2\pi}e^{i(2n-k)x}e^{-ikx}\,\mathrm{d}x\\ &=\frac{2\pi}{4^n}\binom{2n}{n}\tag{1a} \end{align} $$ In this answer, it is shown geometrically that for $0\le x\lt\frac\pi2$, we have $\sin(x)\le x\le\tan(x)$.

Since $\sin^2(x)\le x^2$, $\cos^2(x)=1-\sin^2(x)\ge1-x^2$. Furthermore, $e^{x^2}\ge1+x^2$. Thus, $$ \cos^2(x)\,e^{x^2}\ge1-x^4\tag{1b} $$ For $0\le x\le\frac\pi2$, we have $\frac{\mathrm{d}}{\mathrm{d}x}\cos^2(x)\,e^{x^2}=2(x-\tan(x))\cos^2(x)\,e^{x^2}\le0$. Therefore, $$ \cos^2(x)\,e^{x^2}\le1\tag{1c} $$ $\text{(1b)}$ and $\text{(1c)}$ give $$ 1-x^4\le\cos^2(x)\,e^{x^2}\le1\tag{1d} $$ Using $[a,b]$ to denote a number between $a$ and $b$, Bernoulli's Inequality, $\text{(1d)}$, and $x^4\le\frac{x^3+x^5}2$ yield $$ \begin{align} \sqrt{n}\int_0^{2\pi}\cos^{2n}(x)\,\mathrm{d}x &=4\sqrt{n}\int_0^{\pi/2}\cos^{2n}(x)\,\mathrm{d}x\\ &=4\sqrt{n}\int_0^{\pi/2}\left[\cos^2(x)\,e^{x^2}\right]^{\normalsize n}e^{-nx^2}\,\mathrm{d}x\\ &=4\sqrt{n}\int_0^{\pi/2}\left[1-nx^4,1\right]e^{-nx^2}\,\mathrm{d}x\\ &=4\int_0^{\pi\sqrt{n}/2}\left[1-\frac{x^4}n,1\right]e^{-x^2}\,\mathrm{d}x\\ &=4\int_0^{\pi\sqrt{n}/2}e^{-x^2}\,\mathrm{d}x-\left[0,\frac3n\right]\tag{1e} \end{align} $$ The Squeeze Theorem applied to $\text{(1e)}$ gives $$ \begin{align} \lim_{n\to\infty}\sqrt{n}\int_0^{2\pi}\cos^{2n}(x)\,\mathrm{d}x &=4\int_0^\infty e^{-x^2}\,\mathrm{d}x\\[6pt] &=2\sqrt\pi\tag{1f} \end{align} $$ Combining $\text{(1a)}$ and $\text{(1f)}$ yields $(1)$.

${\large\square}$

---

Theorem $\bf{2}$: $$ \log\left(\frac{n!\,e^n}{\sqrt{n}\,n^n}\right)=1-\sum_{k=1}^{n-1}\int_0^1\frac{\left(x-\frac12\right)^2}{k^2+k+x(1-x)}\,\mathrm{d}x\tag2 $$

Theorems $3$ and $4$ follow fairly simply from Theorems $1$ and $2$.

Proof: ( Euler-Maclaurin applied to $\log(n!)$ ) $$ \begin{align} \log(n!) &=\sum_{k=1}^n\log(k)\tag{2a}\\ &=\int_1^{n^+}\log(x)\,\mathrm{d}\lfloor x\rfloor\tag{2b}\\ &=\int_1^n\log(x)\,\mathrm{d}x-\int_1^{n^+}\log(x)\,\mathrm{d}\!\left(\{x\}-\tfrac12\right)\tag{2c}\\ &=n\log(n)-n+1+\tfrac12\log(n)+\int_1^n\frac{\{x\}-\frac12}x\,\mathrm{d}x\tag{2d}\\ &=\log\left(\sqrt{n}\,\frac{n^n}{e^n}\right)+1+\sum_{k=1}^{n-1}\int_0^1\frac{x-\frac12}{k+x}\,\mathrm{d}x\tag{2e}\\ &=\log\left(\sqrt{n}\,\frac{n^n}{e^n}\right)+1-\sum_{k=1}^{n-1}\int_0^1\frac{\left(x-\frac12\right)^2}{\left(k+\frac12\right)\left(k+x\vphantom{\frac12}\right)}\,\mathrm{d}x\tag{2f}\\ &=\log\left(\sqrt{n}\,\frac{n^n}{e^n}\right)+1-\sum_{k=1}^{n-1}\int_0^1\frac{\left(x-\frac12\right)^2}{k^2+k+x(1-x)}\,\mathrm{d}x\tag{2g}\\ \end{align} $$ Explanation:
$\text{(2a)}$: write the log of a product as the sum of logs
$\text{(2b)}$: apply the Riemann-Stieltjes integral
$\text{(2c)}$: $\lfloor x\rfloor=x-\{x\}$
$\text{(2d)}$: evaluate the left integral and integrate by parts on the right
$\text{(2e)}$: collect terms and write the integral with $\{x\}$ as a sum of integrals on $[0,1]$
$\text{(2f)}$: subtract $\int_0^1\frac{x-1/2}{k+1/2}\,\mathrm{d}x=0$ from each term in the sum
$\text{(2g)}$: substitute $x\mapsto1-x$ and average

${\large\square}$

---

Theorem $\bf{3}$: (Stirling's Asymptotic Approximation) $$ \lim_{n\to\infty}\frac{n!\,e^n}{\sqrt{2\pi n}\,n^n}=1\tag3 $$

Proof: Define $a_n=\frac{n!\,e^n}{\sqrt{n}\,n^n}$. $(2)$ says that $a_n$ is decreasing. Because $\frac1{k^2+k+x(1-x)}\le\frac1{k(k+1)}$ and $\int_0^1\left(x-\frac12\right)^2\,\mathrm{d}x=\frac1{12}$, $(2)$ also says that $$ \begin{align} \log(a_n) &\ge1-\frac1{12}\sum_{k=1}^\infty\frac1{k(k+1)}\\ &=\frac{11}{12}\tag{3a} \end{align} $$ Therefore, $\lim\limits_{n\to\infty}a_n$ exists and is positive. Finally, $$ \begin{align} 1 &=\lim_{n\to\infty}\frac{\sqrt{\pi n}}{4^n}\binom{2n}{n}\tag{3b}\\[3pt] &=\lim_{n\to\infty}\frac{\sqrt{\pi n}}{4^n}\frac{(2n)!}{n!^2}\tag{3c}\\ &=\lim_{n\to\infty}\frac{\sqrt{\pi n}}{4^n}\frac{a_{2n}\sqrt{2n}\frac{(2n)^{2n}}{e^{2n}}}{\left(a_n\sqrt{n}\frac{n^n}{e^n}\right)^2}\tag{3d}\\ &=\sqrt{2\pi}\lim_{n\to\infty}\frac{a_{2n}}{a_n^2}\tag{3e}\\[6pt] &=\sqrt{2\pi}\big/\lim_{n\to\infty}a_n\tag{3f} \end{align} $$ Explanation:
$\text{(3b)}$: Theorem $1$
$\text{(3c)}$: write binomial coefficients using factorials
$\text{(3d)}$: apply the definition of $a_n$
$\text{(3e)}$: cancel equal terms
$\text{(3f)}$: since $\lim\limits_{n\to\infty}a_n$ exists and is not $0$, it is the reciprocal of $\lim\limits_{n\to\infty}\frac{a_{2n}}{a_n^2}$

$(3)$ is the reciprocal of $\text{(3f)}$.

${\large\square}$

---

Theorem $\bf{4}$: (Bounds on Stirling's Formula) $$ 1+\frac1{12n+1}\le\frac{n!\,e^n}{\sqrt{2\pi n}\,n^n}\le1+\frac1{12n-1}\tag4 $$

Proof: Combining Theorem $2$ and Theorem $3$ yields $$ \log\left(\frac{n!\,e^n}{\sqrt{2\pi n}\,n^n}\right) =\sum_{k=n}^\infty\int_0^1\frac{\left(x-\frac12\right)^2}{k^2+k+x(1-x)}\,\mathrm{d}x\tag{4a} $$ Since $\int_0^1\left(x-\frac12\right)^2\,\mathrm{d}x=\frac1{12}$, for $0\le x\le1$ and $k\ge1$, $$ \frac1{12}\frac1{\left(k+\frac1{12}\right)\left(k+\frac{13}{12}\right)}\le\int_0^1\frac{\left(x-\frac12\right)^2}{k^2+k+x(1-x)}\,\mathrm{d}x\le\frac1{12}\frac1{k(k+1)}\tag{4b} $$ Plugging $\text{(4b)}$ into $\text{(4a)}$ gives $$ \frac1{12n+1}\le\log\left(\frac{n!\,e^n}{\sqrt{2\pi n}\,n^n}\right)\le\frac1{12n}\tag{4c} $$ Noting that $1+x\le e^x\le\frac1{1-x}$ leads to $(4)$.

${\large\square}$

Answer 11

The short and intuitive proof of Stirling formula I was taught long time ago relies on steepest descent/saddle point evaluation of Gamma integral. Let us expand $f(x)=a\log x-x$ in Taylor series around its maximum $x_m$: $0=f'(x)=a/x-1\Rightarrow x_m=a$, $f''(x)=-a/x^2$ hence for $a\rightarrow\infty$, $\Gamma(a+1)=\int_0^\infty e^{f(x)}dx\rightarrow e^{a\log{a}-a}\int_{-\infty}^{+\infty}e^{-(x-x_m)^2/(2a)}dx = e^{a\log{a}-a}\sqrt{2\pi a}$. Exactness of the proof depends on demonstration that contributions to integral away from $x_m$ become insignificant. For me it suffices that $\sigma/x_m\equiv\sqrt{a}/a\rightarrow 0$ for $a\rightarrow\infty$ but more elaborate discussion is possible. This relies on Gauss integral, proven e.g. as $\left[\int_{-\infty}^{+\infty}e^{-x^2}dx\right]^2=\int_{-\infty}^{+\infty}e^{-x^2-y^2}dx\,dy=\int_0^\infty e^{-r^2}2\pi r\,dr=\pi$

Answer 12

Using the Leibniz rule $\frac{d}{dt}\Bigg[\int_{a(t)}^{b(t)}f(x,t)dx\Bigg]=\int_{a(t)}^{b(t)}\frac{\partial}{\partial t}f(x,t)dx+\bigg[f(b(t),t)\frac{d}{dt}b(t)-f(a(t),t)\frac{d}{dt}a(t)\bigg]$

$$ \int_0^\infty e^{-x}dx=1\\ \text{Let }x=tu\implies dx=t.du,t>0\\ \int_0^\infty te^{-tu}du=1\implies \int_0^\infty e^{-tu}du=\frac{1}{t}\\ \frac{d}{dt}\int_0^\infty e^{-tu}du=\int_0^\infty \frac{\partial}{\partial t}e^{-tu}du=\frac{-1}{t^2}\\ \int_0^\infty ue^{-tu}du=\frac{1}{t^2}\\ \int_0^\infty u^2e^{-tu}du=\frac{2}{t^3}\\ \int_0^\infty u^3e^{-tu}du=\frac{6}{t^4}\\ \int_0^\infty u^4e^{-tu}du=\frac{24}{t^5}\\ \int_0^\infty u^ne^{-tu}du=\frac{n!}{t^{n+1}}\\ \implies \int_0^\infty \bigg[\frac{x}{t}\bigg]^ne^{-x}\frac{dx}{t}=\frac{n!}{t^{n+1}}\\ \implies \boxed{\int_0^\infty x^ne^{-x}dx=n!=\Gamma(n+1)} $$ Dividing by $\bigg[\frac{n}{e}\bigg]^n$ $$ \bigg[\frac{e}{n}\bigg]^n\Gamma(n+1)=\int_0^\infty \left[\frac{x}{n}\right]^ne^{-(x-n)}dx\\ \text{Let } s=x-n,\\ \bigg[\frac{e}{n}\bigg]^n\Gamma(n+1)=\int_{-n}^\infty \bigg[1+\frac{s}{n}\bigg]^ne^{-s}ds=\int_{-n}^\infty f(s)ds $$ where $f(s)=\bigg[1+\frac{s}{n}\bigg]^ne^{-s}$ $$ \ln(f(s))=n\ln\bigg[1+\frac{s}{n}\bigg]-s=n\bigg[\frac{s}{n}-\frac{1}{2}\bigg(\frac{s}{n}\bigg)^2+\frac{1}{3}\bigg(\frac{s}{n}\bigg)^3-\cdots\bigg]-s\\ =\frac{-s^2}{2n}+\frac{s^3}{3n^2}-\cdots $$ For large $n$, $$ \ln(f(s))\approx\frac{-s^2}{2n} $$ $$ \bigg(\frac{e}{n}\bigg)^n\Gamma(n+1)=\int_{-\infty}^\infty e^{-s^2/2n}ds=\sqrt{2\pi n}, \text{ since }\int_{-\infty}^\infty e^{-at^2}dt=\sqrt{\frac{\pi}{a}}\\ \Gamma(n+1)=\bigg(\frac{n}{e}\bigg)^n\sqrt{2\pi n}\\ \implies \boxed{n!=\Gamma(n+1)\approx \bigg(\frac{n}{e}\bigg)^n\sqrt{2\pi n}} $$