Why is Euler's Gamma function the "best" extension of the factorial function to the reals?

There are lots (an infinitude) of smooth functions that coincide with $f(n)=n!$ on the integers. Is there a simple reason why Euler's Gamma function $\Gamma (z) = \int_0^\infty t^{z-1} e^{-t} dt$ is the "best"? In particular, I'm looking for reasons that I can explain to first-year calculus students.

The Bohr–Mollerup theorem shows that the gamma function is the only function that satisfies the properties

• $f(1)=1$;
• $f(x+1)=xf(x)$ for every $x\geq 0$;
• $\log f$ is a convex function.

The condition of log-convexity is particularly important when one wants to prove various inequalities for the gamma function.

By the way, the gamma function is not the only meromorphic function satisfying $$f(z+1)=z f(z),\qquad f(1)=1,$$ with no zeroes and no poles other than the points $z=-n$, $n=0,1,2\dots$. There is a whole family of such functions, which, in general, have the form $$f(z)=\exp{(-g(z))}\frac{1}{z\prod\limits_{m=1}^{\infty} \left(1+\frac{z}{m}\right)e^{-z/m}},$$ where $g(z)$ is an entire function such that $$g(z+1)-g(z)=\gamma+2k\pi i,\quad k\in\mathbb Z, $$ ($\gamma$ is Euler's constant). The gamma function corresponds to the simplest choice $g(z)=\gamma z$.

Edit: corrected index in the product.

Actually there are other (less-frequently) used extensions to the factorial, with different properties from the gamma function which may be desirable in some contexts.

Euler's Gamma Function

Hadamard's Gamma function

Luschny's factorial function

See here for more information.

For whatever reason, Nature (by which I mean integrals) seems to prefer the Gamma function as the "correct" substitute for the factorial in various integrals, which seems to come more or less from its integral definition. For example, for non-negative integers $a, b$, it's not hard to show (and there's a really nice probabilistic argument) that

$\displaystyle \int_{0}^1 t^a (1 - t)^b \, dt = \frac{a! b!}{(a+b+1)!}.$

For (non-negative?) real values of $a$ and $b$ the correct generalization is

$\displaystyle \int_0^1 t^a (1 - t)^b \, dt = \frac{\Gamma(a+1) \Gamma(b+1)}{\Gamma(a+b+2)}.$

And, of course, integrals are important, so the Gamma function must also be important. For example, the Gamma function appears in the general formula for the volume of an n-sphere. But the result that, for me, really forces us to take the Gamma function seriously is its appearance in the functional equation for the Riemann zeta function.

Wielandt's theorem says that the gamma-function is the only function $f$ that satisfies the properties:

• $f(1)=1$
• $f(z+1)=zf(z)$ for all $z>0$
• $f(z)$ is analytic for $\operatorname{Re}z>0$
• $f(z)$ is bounded for $1\leq \operatorname{Re}z\leq 2$

(See also the related MathOverflow thread Importance of Log Convexity of the Gamma Function, where I learned about the above theorem.)