Собственно, интеграл в правой части называется гамма-функцией — и в него можно подставлять не только целые значения показателя у x:
Г(a) := \int_0^{\infty} x^{a-1} e^{-x} dx;
так что формула выше утверждает, что n!=Г(n+1).

Интегрируя по частям, легко проверить, что Г(a+1)=a Г(a), и по индукции проверить, что действительно
Г(n+1)=n!.

Так вот — а как бы нам оценить, чему равен этот интеграл?

Сначала посмотрим: а где подынтегральная функция f(x)=x^n e^{-x} принимает максимальное значение?

Запишем её как exp(g(x)), где
g(x)= ln f(x) = (n ln x - x);

тогда мы хотим максимизировать g(x). Её производная равна g'(x)=n/x -1, значит, максимальное значение принимается в точке x_0=n.

И равно оно f(n)=(n/e)^n — то есть мы опять поймали экспоненциальную часть приближения.

Вынесем его за интеграл; останется интеграл от
exp(g(x)-g(x_0)).

Вот как выглядит график подынтегральной функции при n=30:

Правда, напоминает гауссовский колокол выше?

На самом деле — он и есть (чем больше n, тем точнее).

Потому что — посмотрим, как ведёт себя функция g(x)-g(n) рядом с точкой x=n, где её значение максимально.

Значение, собственно, у неё там обращается в 0 — мы его вычли.

Производной тоже нет: это точка максимума.

Значит, первый нетривиальный член ряда Тейлора квадратичный;
g''(n) = -1/n,
поэтому

А тогда

На самом деле, приближение выше хорошо работает не только при маленьких y, но при y порядка корня из n — а потом подынтегральная функция уже становится слишком маленькой.

Поэтому гамма-функция Г(n) примерно равна произведению (n/e)^n на интеграл от e^{-y^2/2n}.