问题起因
阶乘 \(n!\) 的增长速度非常快。\(20!\) 不能存储在典型的 int
变量中,\(200!\) 就连双精度浮点变量也不能近似。处理阶乘的对数会是更方便的选择。
那么,该如何在不计算阶乘结果的前提下,计算阶乘的对数?
斯特林公式
斯特林公式(Stirling's approximation)是一条用来取 \(n!\) 的近似值的数学公式。斯特林公式能够将求解阶乘的复杂度降低到对数级。
\[n!\approx \sqrt{2\pi n}(\frac{n}{e})^n \]其误差约为 \(\frac{1}{12n}\)。换句话说,\(n\) 越大误差越小。
借助斯特林公式,可以方便地计算阶乘的对数。
\[\ln n! \approx (n-\frac{1}{2})\ln n -n +\frac{1}{2}\ln 2\pi \]提高精度
Gergő Nemes在2007年提出了一个近似公式:
\[\ln n!\approx \frac{1}{2}[\ln (2\pi)-\ln n]+n\left[\ln\left(n+\frac{1}{12n-\frac{1}{10n}}\right)-1\right] \]当 \(n\) 大于 8 时,近似值可精确到小数点后 8 位。
Python 代码与一个奇特技巧
查阅 radarsimpy 的代码时,发现有一个取巧的小想法。
既然 \(n\) 越大时估算越准确,那么为何不计算 \(n+x\) 的结果,然后把 \(x\) 给去掉?
def log_factorial(n):
n = n + 9.0 # 让 n 变大
product = 1
for i in range(1, 9):
product *= n - i
return (
1 / 2 * (np.log(2 * np.pi) - np.log(n))
+ n * (np.log(n + 1 / (12 * n - (1 / 10 / n))) - 1)
- np.log(product) # 去掉多算的阶
)
参考来源
- John,“How to compute log factorial”,https://www.johndcook.com/blog/2010/08/16/how-to-compute-log-factorial/
- https://zh.wikipedia.org/wiki/史特靈公式