The factorial recursion \(n!=n\cdot(n-1)!\) is exact but tells you almost nothing about how fast \(n!\) grows — computing it for large \(n\) just produces an opaque number. What is missing is a closed-form expression that tracks \(n!\) closely as \(n\to\infty\). Building one from scratch, starting from the plain definition \(n!=1\cdot2\cdot3\cdots n\), turns out to explain every piece of the final formula, including the otherwise mysterious appearance of \(\pi\).
From a Product to a Sum
Products are awkward to approximate directly, because a small relative error in one factor multiplies into every factor after it. Taking a logarithm turns the product into a sum, where errors only add:
\[\ln(n!)=\ln1+\ln2+\cdots+\ln n=\sum_{k=1}^{n}\ln k.\]
Plotting the terms \(\ln1,\ln2,\ldots,\ln n\) as the heights of unit-width bars, this sum is exactly the total area of those bars. Since \(\ln x\) is increasing, that bar area is close to the area under the smooth curve \(y=\ln x\) itself, suggesting the sum can be approximated by the corresponding integral:
\[\sum_{k=1}^{n}\ln k\approx\int\limits_1^n\ln x\,dx.\]
Integration by parts with \(u=\ln x\), \(dv=dx\) gives \(\int\ln x\,dx=x\ln x-x\), so
\[\int\limits_1^n\ln x\,dx=\bigl[x\ln x-x\bigr]_1^n=n\ln n-n+1.\]
Exponentiating \(\ln(n!)\approx n\ln n-n+1\) already produces the dominant part of the final formula:
\[n!\approx e\left(\frac{n}{e}\right)^n.\]
This pins down where the \((n/e)^n\) term and the base of the natural logarithm come from: they fall directly out of the antiderivative of \(\ln x\), long before any mention of Gaussian integrals or Wallis products. What is still missing is everything the crude "sum equals integral" step glossed over.
A Half-Power Correction
Approximating the sum by the integral treats every bar as if its area matched the curve underneath it exactly, which overstates the area on the rising parts of \(\ln x\) and understates it near the bar tops. The standard fix is the trapezoidal rule: instead of the plain integral, average in half of each endpoint's function value,
\[ \sum_{k=1}^{n}\ln k \approx \int\limits_1^n\ln x\,dx+\frac{\ln1+\ln n}{2} =n\ln n-n+1+\frac{\ln n}{2}, \]
using \(\ln1=0\). The trapezoidal correction is not itself exact either — matching curvature, not just slope, requires yet more terms — so the constant contribution from this whole family of corrections is best left as a single unknown constant \(C\) for now rather than trusted as the literal \(+1\) above:
\[\ln(n!)=n\ln n-n+\tfrac12\ln n+C+o(1),\qquad n\to\infty.\]
Exponentiating introduces the half-integer power of \(n\) as \(e^{\frac12\ln n}=\sqrt n\):
\[n!\sim e^C\sqrt n\left(\frac ne\right)^n=A\sqrt n\left(\frac ne\right)^n,\qquad A:=e^C.\]
So \(\sqrt n\) is not a separate insight bolted onto the formula — it is exactly the boundary term the trapezoidal rule adds on top of the plain rectangle-integral approximation. What remains is to pin down the one constant, \(A\), that the local sum-versus-integral comparison cannot see.
Higher-Order Corrections
Repeating the "sum equals integral plus boundary corrections" idea to higher order — comparing a sum to an integral using successively higher derivatives of the summand at the endpoints — is systematized by the Euler–Maclaurin formula (Leonhard Euler and Colin Maclaurin, independently, around 1735). Applied to \(f(x)=\ln x\) on \([1,n]\), it extends the expansion above with an asymptotic series in descending powers of \(n\):
\[ \ln(n!)=\left(n+\tfrac12\right)\ln n-n+C+\frac{1}{12n}-\frac{1}{360n^3}+\frac{1}{1260n^5}-\cdots. \]
Exponentiating turns the correction terms into a multiplicative factor:
\[ n!\sim A\sqrt n\left(\frac ne\right)^n\exp\left(\frac{1}{12n}-\frac{1}{360n^3}+\cdots\right). \]
Expanding that exponential itself as a Taylor series in \(1/n\) reproduces the familiar polynomial-looking correction factor
\[ n!\sim A\sqrt n\left(\frac ne\right)^n\left(1+\frac{1}{12n}+\frac{1}{288n^2}-\frac{139}{51840n^3}+\cdots\right). \]
Every coefficient here is generated mechanically from the Euler–Maclaurin expansion of \(\ln x\). The one thing this entire local analysis cannot determine is the constant \(A\) itself, since a comparison of a sum to an integral near \(x=n\) has no way of seeing back to the fixed starting point \(x=1\). Finding \(A\) needs an independent identity where the same unknown constant shows up twice and can be solved for.
Pinning Down the Constant with Wallis' Product
The central binomial coefficient \(\binom{2n}{n}=\dfrac{(2n)!}{(n!)^2}\) supplies exactly such an identity, because its large-\(n\) behavior is already known independently through Wallis' product (John Wallis, 1656), a classical infinite product for \(\pi\) built from ratios of even and odd integers. Its consequence for the central binomial coefficient is the asymptotic estimate
\[\binom{2n}{n}\sim\frac{4^n}{\sqrt{\pi n}},\qquad n\to\infty.\]
Substituting the still-unknown-constant form \(n!\sim A\sqrt n\,(n/e)^n\) into both \((2n)!\) and \((n!)^2\) gives
\[ (2n)!\sim A\sqrt{2n}\left(\frac{2n}{e}\right)^{2n}, \qquad (n!)^2\sim A^2n\left(\frac ne\right)^{2n}. \]
Dividing the first by the second, the \((n/e)^{2n}\) factors combine into a clean power of \(4\):
\[ \binom{2n}{n}=\frac{(2n)!}{(n!)^2} \sim \frac{A\sqrt{2n}}{A^2n}\cdot\frac{(2n)^{2n}}{n^{2n}} =\frac{\sqrt2}{A}\cdot\frac{4^n}{\sqrt n}. \]
Matching this against the Wallis-product estimate \(4^n/\sqrt{\pi n}\) forces the prefactors to agree:
\[\frac{\sqrt2}{A}=\frac{1}{\sqrt\pi}\quad\Longrightarrow\quad A=\sqrt{2\pi}.\]
With the constant fixed, the sum-versus-integral derivation is complete:
\[\boxed{n!\sim\sqrt{2\pi n}\left(\frac ne\right)^n}\qquad(n\to\infty).\]
A Second Route: Laplace's Method on the Defining Integral
The Gamma function represents the factorial not as a product but as the integral \(n!=\int_0^{\infty}x^ne^{-x}\,dx=\int_0^{\infty}e^{n\ln x-x}\,dx\). Reaching the same asymptotic formula from there needs none of the discrete sum machinery above — only Laplace's method, the general technique of approximating an integral by expanding its exponent around the single point where it is largest.
Substituting \(x=n(1+t)\), so \(dx=n\,dt\), turns the exponent into
\[n\ln x-x=n\ln n-n+n\bigl[\ln(1+t)-t\bigr],\]
and the integral into
\[n!=n^ne^{-n}\int\limits_{-1}^{\infty}e^{\,n[\ln(1+t)-t]}\,n\,dt.\]
Near \(t=0\), Taylor expansion gives \(\ln(1+t)-t=-\tfrac12t^2+O(t^3)\), and for large \(n\) the factor \(e^{n[\ln(1+t)-t]}\) is negligible everywhere except in a shrinking neighborhood of \(t=0\), so replacing the bracket by its leading quadratic term and extending the limits of integration to \(\pm\infty\) changes the integral by an amount that vanishes relative to the whole:
\[ n!\approx n^{n+1}e^{-n}\int\limits_{-\infty}^{\infty}e^{-nt^2/2}\,dt =n^{n+1}e^{-n}\sqrt{\frac{2\pi}{n}} =\sqrt{2\pi n}\left(\frac ne\right)^n, \]
using the same Gaussian integral that defines \(\int_{-\infty}^{\infty}e^{-x^2}\,dx=\sqrt\pi\), rescaled by \(\sqrt{n/2}\). That this second, completely independent route — steepest descent on a continuous integral instead of matching a discrete binomial identity — lands on exactly the same constant \(\sqrt{2\pi}\) is a good consistency check on both derivations: \(\pi\) enters here through the Gaussian bell curve, and in the previous section through Wallis' product, yet the two unrelated appearances of \(\pi\) agree to the letter.
The Refined Approximation and a Worked Example
Already at \(n=10\), the leading term alone gives
\[ \sqrt{2\pi\cdot10}\left(\frac{10}{e}\right)^{10}\approx3{,}598{,}696 \]
against the exact value \(10!=3{,}628{,}800\) — a relative error under \(1\%\). Including the first correction term from the Euler–Maclaurin expansion,
\[ n!\sim\sqrt{2\pi n}\left(\frac ne\right)^n\left(1+\frac{1}{12n}\right), \]
brings the same \(n=10\) estimate to within \(0.003\%\) of the exact factorial. Unlike a Taylor series, the full expansion in descending powers of \(n\) diverges for every fixed \(n\); it is an asymptotic series, useful precisely because truncating it after finitely many terms keeps improving the approximation up to a point, before the omitted terms start growing again.
Why \(\sqrt{2\pi}\) Is the Only Constant That Works
Both derivations above fixed the leading constant at \(\sqrt{2\pi}\), but it is worth checking that no other constant could have worked just as well within the same approximation shape \(n!\approx C\sqrt n\,(n/e)^n\). Dividing the exact factorial by this candidate approximation and taking the limit,
\[ \lim_{n\to\infty}\frac{n!}{C\sqrt n\,(n/e)^n} =\frac{1}{C}\lim_{n\to\infty}\frac{n!}{\sqrt n\,(n/e)^n} =\frac{\sqrt{2\pi}}{C}, \]
using the boxed result from above for the remaining limit. This ratio tends to \(1\) — meaning the relative error of the approximation vanishes as \(n\to\infty\) — if and only if \(C=\sqrt{2\pi}\). Any other choice leaves a fixed, nonzero relative error no matter how large \(n\) gets: taking \(C=1\), for instance, gives
\[\frac{n!}{\sqrt n\,(n/e)^n}\to\sqrt{2\pi}\approx2.5066,\]
so that approximation would always overshoot the true factorial by a factor of roughly \(2.51\), for every \(n\), however large. The plot below tracks the exact ratio \(n!/\bigl(\sqrt{2\pi n}(n/e)^n\bigr)\) for small \(n\): it starts above \(1\) and decreases monotonically toward the dashed line at \(y=1\), visibly approaching — but for any finite \(n\) never reaching — the constant found above.