Where we left the gamma function
In the previous guide you built the gamma function from the Euler integral Gamma(z) = integral from 0 to infinity of t^(z-1) e^(-t) dt, which converges for Re z > 0, and you saw it satisfies the functional equation Gamma(z + 1) = z Gamma(z). That recursion is the workhorse: it interpolates the factorial (Gamma(n + 1) = n!) and, run backwards as Gamma(z) = Gamma(z + 1)/z, it pushes Gamma to the whole plane as a meromorphic function, with simple poles at z = 0, -1, -2, ... and no zeros anywhere.
Two facts from that guide are the seeds of everything here. First, Gamma is zero-free — its reciprocal 1/Gamma(z) is an entire function, one that vanishes exactly at the non-positive integers. Second, by Weierstrass factorization an entire function is built from its zeros, so 1/Gamma must be expressible as a product over the points 0, -1, -2, .... Chase that thought and you get the Weierstrass product for the gamma function — the bridge to the reflection formula.
Written out, that product is 1/Gamma(z) = z e^(gamma z) times the product over n >= 1 of (1 + z/n) e^(-z/n), where gamma = 0.5772... is the Euler-Mascheroni constant. The bare factor z handles the zero at the origin; each factor (1 + z/n) contributes the zero at z = -n; and the e^(-z/n) is precisely the Weierstrass correction that makes the slowly-decaying product converge. The constant gamma is forced on us — it is exactly what the e^(-z/n) corrections leave behind, and it shows up the moment you compare growth rates of the two sides.
The reflection formula: gamma meets the sine
Here is the magic trick. We have a product for 1/Gamma(z) and the same product for 1/Gamma(-z). Multiply the two. The bare factors give z times (-z) = -z^2, and each pair of corrected factors combines: (1 + z/n)(1 - z/n) = 1 - z^2/n^2, while the exponentials e^(-z/n) and e^(+z/n) cancel exactly. What survives is -z^2 times the product over n of (1 - z^2/n^2) — and that product is something you met two rungs ago.
It is the sine product: sin(pi z) = pi z times the product over n of (1 - z^2/n^2). So the product over n of (1 - z^2/n^2) equals sin(pi z)/(pi z), and our computation becomes 1/(Gamma(z) Gamma(-z)) = -z^2 times sin(pi z)/(pi z) = -z sin(pi z)/pi. Now use the functional equation once: Gamma(1 - z) = -z Gamma(-z), so Gamma(z) Gamma(1 - z) = -z Gamma(z) Gamma(-z). Substituting gives the prize.
Gamma(z) * Gamma(1 - z) = pi / sin(pi z) z = 1/2 : Gamma(1/2)^2 = pi / sin(pi/2) = pi => Gamma(1/2) = sqrt(pi)
This is the reflection formula: Gamma(z) Gamma(1 - z) = pi/sin(pi z). Read it as a conservation law. The poles of Gamma(z) at z = 0, -1, -2, ... and the poles of Gamma(1 - z) at z = 1, 2, 3, ... together account for poles at every integer — exactly where pi/sin(pi z) blows up. And since the right side is never zero (sin(pi z) is finite), the left side never vanishes, which re-proves that Gamma has no zeros. One identity, and the entire zero-and-pole bookkeeping of gamma snaps into place.
A few dividends, including the duplication formula
The reflection formula pays out immediately. Setting z = 1/2 gives Gamma(1/2)^2 = pi, hence Gamma(1/2) = sqrt(pi) — the source of all those sqrt(pi)'s in Gaussian integrals. More generally it lets you trade a value of Gamma in the right half plane for one in the left, which is how you compute Gamma at, say, z = -1/2 + i without ever leaving the convergent integral region. It is a genuine reflection across the line Re z = 1/2, the same line that will become the star of the zeta story.
A cousin of the reflection formula is the duplication formula: Gamma(z) Gamma(z + 1/2) = 2^(1 - 2z) sqrt(pi) Gamma(2z). Where reflection relates z to 1 - z, duplication relates z and z + 1/2 to 2z — it tells you how gamma behaves when you double the argument. Both are special cases of the Gauss multiplication formula, and both fall out of the same Weierstrass product by pairing factors cleverly. The lesson is that the product representation, not the integral, is the right tool for proving multiplicative identities.
Stirling: how fast does n! grow?
Now the second secret. The functional equation tells you Gamma jumps by a factor of z at each step, but it does not tell you the size of Gamma(z) for large z. Stirling's formula answers that with stunning precision: for large positive x, Gamma(x + 1) = x! is approximately sqrt(2 pi x) times (x/e)^x. So n! is not merely 'big' — it grows like (n/e)^n times a slowly-varying sqrt-of-n correction. For n = 10 this gives about 3.60 million against the true 3.6288 million, already under 1% off, and the relative error shrinks as n grows.
Where does sqrt(2 pi x) (x/e)^x come from? Look at the integral Gamma(x + 1) = integral from 0 to infinity of t^x e^(-t) dt and write the integrand as e^(x log t - t). For large x the exponent x log t - t has a sharp maximum where its derivative x/t - 1 vanishes, namely at t = x. Everything in the integral piles up near t = x; away from there the integrand is exponentially smaller. This is the central idea of Laplace's method: when an integral is dominated by a single peak, you may replace the exponent by its quadratic Taylor approximation at the peak and integrate the resulting Gaussian.
- Write Gamma(x + 1) = integral of e^(phi(t)) dt with phi(t) = x log t - t, and find the peak: phi'(t) = x/t - 1 = 0 gives t = x.
- Evaluate at the peak: phi(x) = x log x - x, so e^(phi(x)) = (x/e)^x — that is the dominant factor, the 'height' of the bump.
- Approximate the bump by a Gaussian using the curvature phi''(x) = -x/x^2 = -1/x; the width of the peak scales like sqrt(x).
- Integrate the Gaussian: integral of e^(-(t - x)^2/(2x)) dt = sqrt(2 pi x); multiply by the height to get Gamma(x + 1) ~ sqrt(2 pi x) (x/e)^x.
Into the complex plane: saddles and the asymptotic series
The real-variable story is honest but limited: it gives the leading term, and only for real x marching to infinity along the positive axis. The functions we actually care about — gamma in the critical strip, the Riemann zeta function — live in the complex plane, and we want Gamma(z) as |z| goes to infinity in a whole sector, not just along the real line. For that the peak becomes a saddle: a point where the complex exponent has zero derivative, but where the real part rises in one direction and falls in another, like a mountain pass.
This is the saddle point method (or method of steepest descent). You deform the contour of integration — legitimate because the integrand is holomorphic, by Cauchy's theorem — so it crosses the saddle along the direction of steepest descent, where the integrand decays as fast as possible and the local picture is again a clean Gaussian. The leading term reproduces sqrt(2 pi z) (z/e)^z, but now valid for complex z in any sector avoiding the negative real axis. The negative axis is excluded for a familiar reason: that is where the poles of Gamma sit, and log z needs a branch cut.
Push the expansion further and you get the full Stirling series: log Gamma(z) = (z - 1/2) log z - z + (1/2) log(2 pi) + 1/(12 z) - 1/(360 z^3) + .... Here is the crucial honesty: this series diverges for every fixed z. It is an asymptotic series, not a convergent one. Its partial sums first approach log Gamma(z), reaching their best accuracy when you stop after roughly the smallest term, and then the terms turn around and grow. The bigger |z| is, the smaller that best error — but you can never make it exactly zero by summing more terms. Asymptotic and convergent are different promises, and confusing them is a classic trap.