By Long Luo
在科研和工程领域中，阶乘( \(\textit{Factorial}\) )有着广泛的应用。在概率论中，阶乘是计算排列( \(\textit{Permutation}\) )和组合( \(\textit{Combination}\) )时不可或缺的；在物理中，计算粒子系统的状态数以及大型系统的统计分布都要用到阶乘；在计算机中，阶乘则用于图论和组合优化问题。
大家都知道“指数爆炸”这个值，因为指数增长是非常迅猛的。其实阶乘增长要远远快于指数增长，如下图 1 所示为不同算法复杂度增长情况。随着 \(n\) 的增大，阶乘的计算复杂度迅速上升，当处理大 \(n\) 时，计算 \(n!\) 会变得极其复杂且耗时。例如 \(5! = 120\) ，而 \(50! \approx 3.04 \times 10^{64}\) ，这已经是一个非常庞大的数字！直接使用递归方法去求解 \(n!\) 的时间复杂度是 \(O(n)\) ，对于较大的 \(n\) 来说很容易栈溢出。在实际应用中，往往需要计算大数的阶乘，即使目前最先进的计算机去处理极大数的阶乘时，也会面临需要巨大的计算及存储资源消耗问题。
图1. 不同算法时间复杂度 Time Complexity Comparison人们迫切需要找到一种可以快速计算阶乘的方法，在 18 世纪初期，苏格兰数学家詹姆斯·斯特林( \(\textit{James Stirling}\) )提出了斯特林公式 。斯特林公式( \(\textit{Stirling's Approximation}\) )提供了一种近似计算阶乘的方法，特别适用于大 \(n\) 的情况，其标准形式为：
\[ n! \approx {\sqrt {2\pi n}} \,\left( {\frac {n}{e}} \right )^n \tag{1.1} \label{1.1} \]
其对数形式为：
\[ \ln (n!) \approx n \ln n - n \]
这个公式最早是亚伯拉罕·棣莫弗( \(\textit{Abraham de Moivre}\) )在研究二项分布时，为了解决大数阶乘问题时发现的，其形式为：
\[ n! \approx C n^{n + \frac {1}{2}}e^{-n} \tag{1.2} \label{1.2} \]
其中 \(C\) 为某个常量值，斯特林证明了公式中 \(C = \sqrt {2 \pi}\) ，于是冠名权就给了斯特林。
斯特林公式可以大大简化阶乘的计算，特别是当 \(n\) 很大时，它提供了一个非常精确的近似值。斯特林公式使得复杂的阶乘计算可以通过较为简单的幂函数、指数函数和根号运算来完成。相比于直接计算阶乘，它极大地减少了计算量，是大数问题中不可或缺的工具。
斯特林公式中集合了圆周率 \(\pi\) 和自然常数 \(e\) ，这 \(2\) 个数学中最重要的常数，十分独特且具有美感。因为 \(e\) 意味着连续增长，而阶乘就是连续自然数的相乘。而出现 \(\pi\) 的时候，就要问自己 “Where is the circle?”，那么阶乘是如何和几何中的圆扯上关系了呢？
关于斯特林公式的证明，常见的证明是对数形式的证明或者利用伽马函数( \(\textit{Gamma Function}\) )来证明，这里将介绍一种从零开始更易理解的推导方式。
从零开始推导斯特林公式假设我们现在不知道斯特林公式，但我们就是要计算一些大数的阶乘，应该怎么办呢？
面对这样的问题，我们自然会问：有没有一种方法能够近似表示大数阶乘的值？如果能找到一个简单但精确的近似表达式，不仅能够帮助我们理解阶乘的增长规律，还可以为复杂公式的推导和计算带来便利。
也就是说我们的目标是找到这样一个函数 \(f(n)\) ，使得：
\[ n! = f(n) + \varepsilon(n) \tag{2.1} \label{2.1} \]
同时满足当 \(n \rightarrow \infty\) ，误差函数 \(\varepsilon(n) \rightarrow 0\) 。
那 \(f(n)\) 会是什么形式呢？我们希望它是初等函数的组合，比如中学时学过的幂函数、对数函数或者三角函数等，这些函数熟悉又容易计算。
显然有：
\[ n! = \prod _{i=1}^{n} i = 1 \cdot 2 \cdot \ldots \cdot (n - 1) \cdot n \le n^n \]
但随着 \(n\) 增大，\(n!\) 明显小于 \(n^n\) ，所以我们需要添加一个平衡函数 \(g(n)\) ，随着 \(n \rightarrow \infty\) ，\(g(n) \rightarrow 0\)。
令 \(f(n) = n^n g(n)\) ，我们的目标就变成了求解函数 \(g(n)\) ，易得：
\[ \begin{aligned} (n+1)! & = (n+1)^{n+1} g(n+1) \\ n! & = n^n g(n) \end{aligned} \]
那么：
\[ \frac {g(n+1)}{g(n)} = \left( 1 + \frac {1}{n} \right)^{-n} \]
等式右边是不是很熟悉，就是自然常数 \(e\) 的表达式：
\[ \lim_{n \rightarrow \infty} \left( 1 + \frac {1}{n} \right)^{-n} = e^{-1} \]
于是我们轻轻松松就得到了 \(g(n) = e^{-n}\) ，综合一下则有：
\[ f(n) = \left ( \frac {n}{e} \right) ^n \tag{2.2} \label{2.2} \]
这个渐进式足够精确吗？上一节，我们得到了函数 \(f(n)\) 的表达式，这个逼近式效果如何呢？我们来简单测试下：
\[ \frac {f(n + 1)}{f(n)} = \left( 1 + \frac {1}{n} \right)^n (n + 1) e^{-1} \tag{3.1} \label{3.1} \]
我们知道 \(e = \lim _{n \to \infty } \left( 1 + \frac {1}{n} \right)^n\) ，对于足够大的 \(n\) ，可以看出上式 \(\eqref{3.1}\) 接近 \((n + 1)\) ，也就是说函数 \(f(n)\) 和阶乘一样满足递归性质的。
但是这个逼近式并不是特别好，为什么呢？
对公式 \(\eqref{3.1}\) 适当变形，则有：
\[ \frac {f(n + 1)/(n + 1)!}{f(n)/n!} = \left( 1 + \frac {1}{n} \right)^n e^{-1} \tag{3.2} \label{3.2} \]
不难得到：
\[ \frac{f(n)}{n!} = e^{-1} \prod_{t = 1}^{n - 1} \left [ \left(1 + \frac{1}{t}\right)^t e^{-1} \right ] \tag{3.3} \label{3.3} \]
容易看出上式 \(\eqref{3.3}\) 右边明显小于 \(1\) ，随着 \(n\) 越大，越趋近于 \(0\) ， 所以 \(n!\) 和 \(n^ne^{-n}\) 并不是同阶的，因此 \(f(n) = n^n e^{-n}\) 的逼近效果并不是很好，我们需要找到更好的渐进式。
更好的渐进式和第一节类似，我们可以再增加一个函数 \(h(n)\) ：
\[ f(n) = n^n e^{-n} h(n) \]
实际上由公式 \(\eqref{3.3}\) 我们知道：
\[ h(n) = e^{-1} \prod_{t = 1}^{n - 1} \left [ \left(1 + \frac{1}{t}\right)^t e^{-1} \right ] \tag{4.1} \label{4.1} \]
如何求 \(h(n)\) 呢？这是一个连乘式，我们知道对数可以将乘法变成加法，两边同取对数，于是：
\[ \ln \left ( h(n) \right ) = -1 + \sum_{t = 1}^{n - 1} \left [t \ln \left(1 + \frac{1}{t} \right) - 1 \right ] \tag{4.2} \label{4.2} \]
由 泰勒级数 ，我们知道
\[ \ln (1+x) = \sum _{n=1}^{\infty } \frac {(-1)^{n+1}}{n}x^{n} = x - \frac {x^{2}}{2} + \frac {x^{3}}{3} - \cdots \]
这里 \(-1 < x \leq 1\) ，代入公式 \(\eqref{4.2}\) 可得：
\[ \ln \left ( h(n) \right ) = -1 + \sum_{t = 1}^{n - 1} \left [-\frac {1}{2t} + \frac {1}{3t^2} - \frac {1}{4t^3} + \cdots \right ] \tag{4.3} \label{4.3} \]
而公式 \(\eqref{4.3}\) 中，显然有当 \(n \to \infty\) 时， \(\sum_{t = 1}^{n - 1} [\dfrac{1}{3t^2} - \dfrac{1}{4t^3} + \cdots]\) 是收敛的。
这里我们直接引用调和级数( \(\textit{Harmonic Series}\) )的结论，则有：
\[ \sum_{t = 1}^{n - 1} - \frac {1}{2t} \approx - \frac {1}{2} \ln (n) = \ln (n^{- \frac {1}{2}}) \]
也就是说当 \(n \to \infty\) 时， \(\dfrac {h(n)}{\sqrt {n}}\) 趋近于一个常数，这里不妨设为 \(C\) 。
至此， \(n! \sim C \sqrt{n} \left( \dfrac{n}{e} \right)^n\) ，这个公式和棣莫弗发现的公式 \(\eqref{1.2}\) 是一致的。
我们现在已经知道这个值就是 \(\sqrt {2\pi}\) ，那么怎么求出这个常数 \(C\) ？怎么和圆扯上关系？
圆在哪里？考虑二项分布式
\[ \binom{2n}{n} = \dfrac{(2n)!}{n!^2} \tag{5.1} \label{5.1} \]
应用我们上一节得到的逼近式结论，则有：
\[ \binom{2n}{n} \sim \frac{4^n}{C \sqrt{n/2}} \tag{5.2} \label{5.2} \]
我们可以构造如下的表达式：
\[ \frac {2}{1} \times \frac {4}{3} \times \frac {6}{5} \times \ldots \times \frac {2n}{2n - 1} = \frac {2^n n!}{\frac {(2n)!}{2^n n!}} = \frac {4^n}{\binom {2n}{n}} \]
即有：
\[ \prod_{n = 1}^{\infty } \frac {2n}{2n - 1} \approx C \sqrt {n/2} \tag{5.3} \label{5.3} \]
同理可得：
\[ \frac {2}{3} \times \frac {4}{5} \times \frac {6}{7} \times \ldots \times \frac {2n}{2n + 1} = \frac {4^n}{\binom{2n}{n}} \times \frac {1}{2n + 1} \]
即：
\[ \prod_{n = 1}^{\infty } \frac {2n}{2n + 1} \approx \frac {C \sqrt{n/2}}{2n} \tag{5.4} \label{5.4} \]
将公式 \(\eqref{5.3}\) 和 \(\eqref{5.4}\) 两边相乘得到：
\[ \frac { 2 \cdot 2 \cdot 4 \cdot 4 \cdot \ldots \cdot 2n \cdot 2n }{ 1 \cdot 3 \cdot 3 \cdot 5 \cdot 5 \cdot \ldots \cdot (2n -1) (2n + 1)} = \frac {C^2}{4} \tag{5.5} \label{5.5} \]
而上式左边就是著名的沃里斯乘积( \(\textit{Wallis Product}\) ) ，
\[ \begin{aligned} \frac {\pi }{2} & = \prod _{n=1}^{\infty }{\frac {4n^2}{4n^2 - 1}} = \prod _{n=1}^{\infty }\left ( {\frac {2n}{2n - 1}} \cdot {\frac {2n}{2n + 1}} \right ) \\ & = {\Big (}{\frac {2}{1}} \cdot {\frac {2}{3}}{\Big )}\cdot {\Big (}{\frac {4}{3}}\cdot {\frac {4}{5}}{\Big )}\cdot {\Big (}{\frac {6}{5}}\cdot {\frac {6}{7}}{\Big )} \cdot \;\cdots \\ \end{aligned} \]
而沃里斯乘积值为 \(\dfrac {\pi}{2}\) ，最终我们得到了 \(C = \sqrt{2 \pi}\) 。
至于沃里斯乘积为什么等于 \(\dfrac{\pi}{2}\) ，和圆有什么关系， 这篇文章就不详细介绍了。3Blue1Brown 专门有个视频 The Wallis product for pi, proved geometrically 讲解这个公式背后的几何直观。
笔者曾经也思考了很久，试图找到这背后的intuition，后来发现当 \(n \rightarrow \infty\) 时， 数轴上每个数都要参与进来，如下图 2 所示：
\[ n! = \prod _{i=1}^{n}i = 1 \cdot 2 \cdot 3 \cdots (n-2) \cdot (n-1) \cdot n \]
图2. 无线延申的数轴 Number Line而intuition则是来自直线就是无限大的圆，圆是无限大的直线，这样就和圆有关了！
图3. 直线可以视为曲率无穷小的圆总结我们从零推导出了斯特林公式，看似简单的阶乘问题居然需要对极限、对数、积分等都深刻数学工具的探索，而它的结果更是让人惊叹：阶乘的增长竟然可以通过一个与圆周率 \(\pi\) 紧密相连的表达式优雅地刻画。这种关联不仅展示了数学内部的深层统一性，也让我们感受到圆的几何之美如何渗透到数量与计算的世界中。斯特林公式犹如一座桥梁，将离散的阶乘与连续的自然规律连接在了一起，令人不禁感叹数学的简洁与深邃。
参考文献Factorial 阶乘Permutation 排列Combination 组合Time Complexity 时间复杂度Travelling salesman problem, TSP 旅行推销员问题James Stirling 詹姆斯·斯特林Abraham De Moivre 亚伯拉罕·棣莫弗Gamma function 伽马函数Stirling’s approximation 斯特林公式Wolfram: Stirling’s ApproximationNatural Number 自然常数 eJohn Napier 纳皮尔Logarithm 对数Taylor Series 泰勒级数Harmonic series 调和级数Wallis product 沃里斯乘积3Blue1Brown: The Wallis product for pi, proved geometricallyRaffi Grinberg:《普林斯顿数学分析读本》Dr. Trefor Bazett: Stirling’s Incredible Approximation Gamma Functions, Gaussians, and Laplace’s MethodFlammable Maths: Factorial’s Asymptotic Expansion - DERIVING STIRLING’S FORMULA!Faculty of Khan: The Stirling Approximation: a 5-minute Derivation!Michel van Biezen: Physics 32.5 Statistical Thermodynamics (7 of 39) Stirling’s Approximation ExplainedPhysics mee: Stirling approximation derivationNumberphile: Big FactorialsMathKiwi: Why is Pi here? | Half factorial without Gamma function