# Stirling's formula: proof?

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$$

-
Keith Conrad has a good explanation of this. google.com/… – jspecter Dec 28 '11 at 16:54
Stirling's formula is a pretty hefty result, so the tools involved are going to go beyond things like routine application of L'Hopital's rule, although I am sure there is a way of doing it that involves L'Hopital's rule as a step. I've just scanned the link posted by jspecter and it looks good and reasonably elementary. Another pretty elementary treatment is Terry Tao's: terrytao.wordpress.com/2010/01/02/…. The two treatments I've taken the time to work through both involved doing residue calculus on the gamma function. – Ben Blum-Smith Dec 28 '11 at 17:10
James: are you interested in proving that C exists or in computing its value? For the former, pretty standard procedures based on an estimation of the ratio of two consecutive terms are enough. – Did Dec 28 '11 at 17:33
There's a proof along this line in The Number $\pi$ by Eymard and Lafon. (A very nice book, btw.) – lhf Dec 31 '11 at 2:12
Dan Romik, Stirling's Approximation for n!: The Ultimate Short Proof?, at JSTOR. – Pierre-Yves Gaillard Jan 3 '12 at 18:13

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

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

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)=\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({t\over 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.

-
Hi Professor Schmuland, can I ask you a quick question regarding Stirling's formula? If I am asked to show that the partial sums of $log(j) = nlogn - n + O(logn)$, and I started off with $\sum log(j) = log(n!) = log(\sqrt{2\pi n}(\frac{n}{e})^n)$ = ... and continued from here, do you think it's ok, and that I am not possibly using the result to prove the result? I am a little concerned about making this first move since on the Wikipedia page of Stirling's formula, the result that I am trying to prove shows up early in the discussion. Thanks for your time, @ByronSchmuland, – User001 Dec 6 '15 at 4:16
+1 Very efficient. -Mark – Dr. MV Mar 25 at 19:30
@Dr.MV Thanks for the kind words. – Byron Schmuland Mar 25 at 19:37
You're welcome. My pleasure. And thank you for this nice post. -Mark – Dr. MV Mar 25 at 19:53

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: $$\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}$$

-
(+1) Nice rewriting of $\sum\limits_{n\ge m}\frac1{n^2}$. That makes for a short derivation of the approximation. My original derivation of the asymptotic series was once much shorter than my current post, but I was nudged to add more details. – robjohn Aug 25 '15 at 14:34
Nice one! ! !(+1) – tired Feb 25 at 1:48

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}

-
Sorry but I do not think I can buy this as a full proof. – Did Dec 31 '11 at 21:32
The main problem is that you apply tools valid for formal power series, to get results on functions coinciding with the sums of power series. Such a move is tempting and there are results which allow it, but in a strict setting only (I recommend pondering the first paragraph of the WP page on formal power series). Imagine for a moment that the recursion is solved by $a_k\sim k^k$. What next? Would you consider that the asymptotics $a_0+a_1/n+a_2/n^2+O(1/n^3)$ .../... – Did Jan 1 '12 at 11:20
.../... is valid? Likewise: what are the dots in (6) and (8)? Some rests of converging series? The next terms of some asymptotic series? Something else? Divergent series are a venerable subject but, here again, some caution is needed. Finally: when you write that Thus, $x$ has a power series in $u$ like $x=0+u+\ldots$, this probably means that, IF $x$ coincides with the sum of such a power series, then the beginning of the power series must be $0+u+$some higher order terms (perfectly correct). But you use it to .../... – Did Jan 1 '12 at 11:20
.../... assert that $x$ IS a power series in $u$ (which is an entirely different matter). // All in all, it might be possible to address these points and to reach a rigorous solution based on this method, but your post does not do that, yet. – Did Jan 1 '12 at 11:20
@Didier: I applied the tools valid for formal power series to get the coefficients of $x$ as a power series in $u$. I was not relying on them to show that $x$ was analytic in $u$; that is handled by Lagrange Inversion. I have now mentioned that explicitly in my answer. The dots in $(6)$ are the next terms in the power series for $x$. The first two sets of dots in $(9)$ are the next terms in the power series, while the last set of dots are the next terms in the asymptotic expansion. I have rearranged and amended my answer; does it now seem more complete? – robjohn Jan 3 '12 at 21:22

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)

-

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^{i\frac{\theta}{\sqrt{n}}}-i\frac{\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.

-
I'd like to have a clarification on the following equation: $$\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)}}$$ Does it use an approximation of the standard normal cumulative distribution function? I can't see which formulation of the central limit theorem is used. Is it a De Moivre-Laplace like formula? – Alessandro Cuttin Mar 25 '13 at 15:13
@ Allessandro Cuttin: No. It is an instance of a Local Central Limit Theorem. These are proved via specral methods, much like the standard proof of the CLT, but the machinery involved is usually heavier. – D. Thomine Aug 11 '13 at 17:56

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.

-

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}$$

-
@James I guess this is what you wanted right? – Pedro Tamaroff Feb 13 '12 at 23:52

Define $\gamma (f_n) := \log f_{n + 1} - \log f_n$. Observe that $f_n \equiv g_n$ if and only if $\gamma (f_n) = \gamma (g_n)$ and $f_1 = g_1$. We have $\gamma (n!) = \log (n + 1)$.

We are going to try to express $n!$ in closed form as $h_n$.

Let

$$h_n = c_n \left ( \frac {n} {e} \right ) ^ n.$$

We have $\gamma (h_n) = \log (n + 1) + \log \frac {c_{n + 1}} {c_n} - \varepsilon_n$ where $c_1 = e$ and $\varepsilon_n \to 0$ is a sequence of positive real numbers. In fact, $\varepsilon_n = 1 - n \log (n + 1)$. Then we have

$$c_n = \frac {n^{\frac {1}{2}} e^n} {\prod_{k = 1}^{n - 1} \left ( 1 + \frac {1} {k} \right )^k} = e^{1 + o (1)} \sqrt n.$$

Hence,

$$\log n! = \left ( \left (n + \frac {1} {2} \right ) \log n - n + 1 \right ) \left ( 1 + o (1) \right ).$$

-
This does not give the right answer since $e\ne\sqrt{2\pi}$. – robjohn Nov 22 '15 at 17:24
I think you should consider that my result is not a strict equality but an asymptotic one. Also, $e = 2.718...$ and $\sqrt {2 \pi} = 2.506...$ is close enough for an asymptotic formula. Anyhow, this is the best I can get using only elementary methods. (Note: Not everyone likes integrals) – Nima Bavari Nov 22 '15 at 18:20
Since the constant is not correct, it is not even an asymptotic approximation. Your result is off by a constant factor, so the ratio of your result and $n!$ does not tend to $1$, which is what an asymptotic approximation does. The statement $$n! = \frac {n^{n + \frac {1}{2}}} {e^{n - 1}} \left ( 1 + o (1) \right )$$ is not true. – robjohn Nov 22 '15 at 18:40
Well, $h_n / n!$ does not tend to 1, but $\log h_n / \log n!$ does, and it was what I meant. What do you suggest I do, remove my answer? – Nima Bavari Nov 22 '15 at 18:46
Perhaps state what you have shown as the final result. That is, $\log(n!)=\log(h_n)(1+o(1))$ – robjohn Nov 22 '15 at 18:49