The factorial \(n!\) is usually introduced by a simple recursion and only for non-negative integers. There is, however, a definite integral that reproduces every factorial exactly, and that same integral still makes sense for numbers that are not integers at all. The route to it starts, somewhat surprisingly, with a completely different landmark integral.
The Gaussian Integral
The Gaussian integral is the statement that the bell curve \(e^{-x^2}\) encloses a finite, exactly known area over the whole real line:
\[\int\limits_{-\infty}^{\infty}e^{-x^2}\,dx=\sqrt{\pi}.\]
As already noted when discussing the Fundamental Theorem of Calculus, \(e^{-x^2}\) has no elementary antiderivative, so this cannot be evaluated by writing down \(F(b)-F(a)\) for a closed-form \(F\). The definite integral over the entire real line has an exact value anyway. Call it \(I\) and square it, writing the second factor with a different name for the integration variable so the product becomes a single double integral over the whole plane:
\[ I^2=\left(\int\limits_{-\infty}^{\infty}e^{-x^2}\,dx\right)\left(\int\limits_{-\infty}^{\infty}e^{-y^2}\,dy\right) =\iint\limits_{\mathbb{R}^2}e^{-(x^2+y^2)}\,dx\,dy. \]
Switching to polar coordinates, \(x=r\cos\theta\), \(y=r\sin\theta\), \(dx\,dy=r\,dr\,d\theta\), turns the exponent \(x^2+y^2\) into the single variable \(r^2\), and the area element gains exactly the factor of \(r\) needed to integrate it in closed form:
\[ I^2 =\int\limits_0^{2\pi}\int\limits_0^{\infty}e^{-r^2}\,r\,dr\,d\theta =2\pi\left[-\tfrac12 e^{-r^2}\right]_0^{\infty} =2\pi\cdot\tfrac12 =\pi. \]
Since the integrand \(e^{-x^2}\) is never negative, \(I\) itself is the positive square root, \(I=\sqrt{\pi}\). This trick — squaring an integral to convert it into a double integral that becomes solvable in polar coordinates — has no analogue for a single-variable substitution and is the reason the Gaussian integral is usually proved this way rather than by finding an antiderivative, which does not exist.
The Factorial
For an integer \(n\ge1\), the ordinary factorial is defined by the finite product and recursion
\[ \begin{array}{rl} n! &= n\cdot(n-1)\cdot(n-2)\cdots1\\ &= n\cdot(n-1)! \end{array} \]
together with the base case \(0!=1\). That base case is usually just declared, since the recursion alone does not produce it. The integral definition introduced next reproduces \(0!=1\) as an actual computation instead of a convention.
From the Factorial Integral to the Gamma Function
Consider the function \(x^2e^{-x}\) on \(x\ge0\). It starts at \(0\), rises to a peak at \(x=2\), and decays back toward \(0\) as the exponential eventually overwhelms the polynomial growth. Shading the area under this curve out to where it has become negligible:
the shaded area works out to \(A=2=2!\). The same construction with \(x^3e^{-x}\), peaking at \(x=3\):
gives \(A=6=3!\). This is not a coincidence of these two particular exponents. For an integer \(n\ge0\), set
\[n!=\int\limits_0^{\infty}x^ne^{-x}\,dx.\]
For \(n\ge1\), integration by parts with \(u=x^n\) and \(dv=e^{-x}\,dx\), so \(du=nx^{n-1}\,dx\) and \(v=-e^{-x}\), gives
\[ \int\limits_0^{\infty}x^ne^{-x}\,dx =\Bigl[-x^ne^{-x}\Bigr]_0^{\infty}+n\int\limits_0^{\infty}x^{n-1}e^{-x}\,dx =n\int\limits_0^{\infty}x^{n-1}e^{-x}\,dx, \]
since \(x^ne^{-x}\to0\) at both ends when \(n\ge1\). This is exactly the recursion \(n!=n\cdot(n-1)!\), now derived from the integral rather than assumed. The case \(n=0\) is evaluated directly rather than by this integration-by-parts step.
Nothing in the integral requires an integer exponent. For a complex number \(z\) with \(\operatorname{Re}(z)>0\), define
\[\boxed{\Gamma(z)=\int\limits_0^{\infty}x^{z-1}e^{-x}\,dx}.\]
Along the positive real axis, \(x^{z-1}=e^{(z-1)\ln x}\) is unambiguous because \(\ln x\) is real. Near \(x=0\), the magnitude behaves like \(x^{\operatorname{Re}(z)-1}\), whose integral converges exactly when \(\operatorname{Re}(z)>0\). As \(x\to\infty\), the exponential factor decays faster than every power, so it guarantees convergence at the other endpoint.
Substituting \(z=n+1\) recovers the factorial integral, hence
\[\Gamma(n+1)=n!\qquad(n\in\mathbb{N}_0).\]
Integration by parts now works for complex \(z\) in the same half-plane and gives the functional equation
\[ \boxed{\Gamma(z+1)=z\,\Gamma(z)}. \]
The base case falls out immediately: \(0!=\int_0^{\infty}e^{-x}\,dx=\bigl[-e^{-x}\bigr]_0^{\infty}=1\), the computation promised above. The half-integer case connects straight back to the Gaussian integral. Substituting \(x=t^2\), \(dx=2t\,dt\) in \(\Gamma(\tfrac12)=\int_0^{\infty}x^{-1/2}e^{-x}\,dx\) gives
\[ \Gamma\!\left(\tfrac12\right) =\int\limits_0^{\infty}t^{-1}e^{-t^2}\cdot2t\,dt =2\int\limits_0^{\infty}e^{-t^2}\,dt =\int\limits_{-\infty}^{\infty}e^{-t^2}\,dt =\sqrt{\pi}, \]
using the evenness of \(e^{-t^2}\) to fold the one-sided integral back into the full Gaussian integral from above. Combined with the recursion, \(0.5!=\Gamma\!\left(\tfrac32\right)=\tfrac12\Gamma\!\left(\tfrac12\right) =\tfrac{\sqrt{\pi}}{2}\).
Properties
With \(\Gamma(z+1)=z\,\Gamma(z)\) and \(\Gamma(1)=0!=1\), unrolling the recursion for a positive integer \(n\) reproduces the ordinary factorial:
\[\Gamma(n)=(n-1)\,\Gamma(n-1)=(n-1)(n-2)\,\Gamma(n-2)=\cdots=(n-1)!,\]
closing the loop back to the familiar factorial and confirming that \(\Gamma\) is exactly the extension it was built to be.
Analytic Continuation
The integral initially defines \(\Gamma(z)\) only for \(\operatorname{Re}(z)>0\), but the functional equation extends it to the left. Choose a positive integer \(m\) such that \(\operatorname{Re}(z+m)>0\). Repeatedly solving the functional equation backward gives
\[ \Gamma(z) =\frac{\Gamma(z+m)}{z(z+1)\cdots(z+m-1)}. \]
The numerator is defined by the convergent integral. The denominator shows exactly where the continuation fails: \(z=0,-1,-2,\ldots\). These points are simple poles. Away from them, the formula agrees with the original integral wherever both definitions apply, so it gives the unique meromorphic continuation of \(\Gamma\) to the complex plane.
The same recursion also pins down exactly how fast \(n!\) grows for large \(n\): a Laplace-method expansion of this very integral around the peak of \(e^{n\ln x-x}\), reusing the Gaussian integral from above, leads to \(n!\sim\sqrt{2\pi n}\,(n/e)^n\), known as Stirling's approximation. The full derivation — including where every piece of that formula comes from and why \(\sqrt{2\pi}\) is the only constant that makes it work — is worked out in Introduction to Stirling's Approximation.