近似计算阶乘的对数

问题起因

阶乘 \(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)  # 去掉多算的阶
    )

不过个人认为还是这样更加实用：

def log_factorial(n):
    if n < 9:
        return np.log(math.factorial(n))
    return (
        1 / 2 * (np.log(2 * np.pi) - np.log(n))
        + n * (np.log(n + 1 / (12 * n - (1 / 10 / n))) - 1)
    )

参考来源

• John，“How to compute log factorial”，https://www.johndcook.com/blog/2010/08/16/how-to-compute-log-factorial/
• https://zh.wikipedia.org/wiki/史特靈公式