Python中的阶乘math.factorial算法

0. 简介

将 $1...n$ 中的所有整数相乘得到一个非负整数 $n$$n$ 的阶乘(阶乘) 并用 $n!$ 表示。但是，我们定义 $0! = 1$。
例如，

4! = 1 \times 2 \times 3 \times 4 = 24

是。

如果我们用 Python 计算，

def factorial_naive(n: int) -> int:
    result = 1
    for i in range(1, n + 1):
        result *= i
    return result

print(factorial_naive(4))  # 24

可以写成它也可以如下计算，如递归函数的介绍中很好的例证。

def factorial_recursive(n: int) -> int:
    if n == 0: return 1
    return n * factorial_recursive(n - 1)

print(factorial_recursive(4))  # 24

如您所见，计算阶乘的代码可以很容易地编写，但是math 模块有一个名为factorial() 的函数，可以更容易地用于计算阶乘。

import math
print(math.factorial(4))  # 24

您可以使用这些方法中的任何一种，但math.factorial() 比factorial_naive() 和factorial_recursive() 更快，并且会产生明显的差异，尤其是当$n$ 很大时。
其实你可以通过分别执行factorial_naive(50000)和math.factorial(50000)来感受不同。

math.factorial() 是用 C 实现的，而 factorial_naive() 是用 Python 编写的，因此您可能期望 math.factorial() 更快。然而，不仅如此，math.factorial() 的算法还可以通过一些技巧来计算比天真的计算更快。

在本文中，我们将了解math.factorial() 如何计算阶乘。

本文使用的源代码可以在下面查看。

https://github.com/python/cpython/blob/3.11/Modules/mathmodule.c

1.方法一：用位移位替换

将整数乘以 $2$ 可以替换为左位移 << 1。由于这成本较低，因此尽可能用位移位代替。
例如，$7!$ 是

\begin{aligned}
7! &= 1 \times 2 \times 3 \times 4 \times 5 \times 6 \times 7 \\
&= 1 \times 2 \times 3 \times 2^2 \times 5 \times (2 \times 3) \times 7 \\
&= 1 \times 3 \times 5 \times 3 \times 7 \times 2^4
\end{aligned}

所以我们知道我们应该首先计算 $1 \times 3 \times 5 \times 3 \times 7 = 315$ 然后位移<< 4。
这将 6 美元的乘法转换为 4 美元的乘法和 1 美元的位移运算。

一般来说

n! = (奇数) \times (2の冪乗)

操作次数可以随着$n$ 的增大而减少。

从现在开始，当$n!$如上表示时，奇数部分将被称为$odd\_part(n)$，$2$的幂的指数部分将被称为$two\_exponent(n)$ . 增加。
例如，当 $n = 6$

6! = (1 \times 3 \times 5 \times 3) \times 2^4

这就是为什么

\begin{aligned}
odd\_part(6) &= 1 \times 3 \times 5 \times 3 = 45\\
two\_exponent(6) &= 4
\end{aligned}

然后可以看到可以计算45 << 4。

从现在开始，我将考虑具体表达 $odd\_part(n)$ 和 $two\_exponent(n)$。

2.制定odd_part

作为一个具体的例子，考虑 $n = 20$ 的情况。
如果我们将整数 $1\dots20$ 按包含的 $2$ 的个数进行分类，则如下所示。

$2$的数量	整数	不包括 $2$ 的数字
$0$	1、3、5、7、9、11、13、15、17、19 美元	1、3、5、7、9、11、13、15、17、19 美元
$1$	2、6、10、14、18 美元	1、3、5、7、9 美元
$2$	4 美元、12 美元、20 美元	1 美元、3 美元、5 美元
$3$	$8$	$1$
$4$	16 美元	$1$

这样看

\begin{aligned}
odd\_part(20) &=(1\times3\times5\times7\times9\times11\times13\times15\times17\times19)\\
&\times(1\times3\times5\times7\times9)\\
&\times(1\times3\times5)\\
&\times(1)\\
&\times(1)\\
\end{aligned}

可以看出
这里正奇数$F(L, U)$ 代表 $L, U$

F(L, U) = \left\{
\begin{aligned}
&[L, U) の奇数の総積 \;\;\;&(L \lt U)\\
&1 & (L \geq U)\\
\end{aligned}
\right. \\

和

\begin{aligned}
odd\_part(20) &= F(1, 21) \times F(1, 11) \times F(1, 7) \times F(1, 3) \times F(1, 3)\\
&= \prod_{i=0}^{4}{F\left(1, \left(\left\lfloor\frac{20}{2^i} \right\rfloor + 1\right) OR\; 1\right)}
\end{aligned}

可以写成其中 $x\; OR\; y$ 表示对整数 $x, y$ 的按位 $or$ 运算。

通常对于正整数 $n$

odd\_part(n) = \prod_{i = 0}^{m}{F\left(1, \left(\left\lfloor\frac{n}{2^i} \right\rfloor + 1\right) OR\; 1\right)}

变成。其中 $m$ 是满足 $\left\lfloor\frac{n}{2^m}\right\rfloor \gt 0$ 的最大整数。

此外，$F(1, U) = F(3, U)$, $F\left(1, \left(\left\lfloor\frac{n}{2^m} \right\rfloor + 考虑到 1 \right) 或\; 1\right) = 1$

odd\_part(n) = \prod_{i = 0}^{m - 1}{F(L_i, U_i)}

然而，

\begin{aligned}
L_i &= 3\\
U_i &= \left(\left\lfloor\frac{n}{2^i} \right\rfloor + 1\right) OR\; 1
\end{aligned}

可以写成

3. 技巧 2：重用预先计算的结果

有了这个，$odd\_part(n)$ 可以用一个公式来表示，但是如果你按原样计算，它不会改变一个幼稚的方法。
因此，通过设计计算顺序来提高效率。
对于正奇数 $a \leq b \leq c$

F(a, c) = F(a, b) \times F(b, c)

使用 $odd\_part(20)$ 变成

\begin{aligned}
odd\_part(20) =& F(3, 21) \times F(3, 11) \times F(3, 7) \times F(3, 3)\\
=& F(3, 3) \times F(3, 7) \times F(7, 11) \times F(11, 21) \times\\
& F(3, 3) \times F(3, 7) \times F(7, 11) \times\\
& F(3, 3) \times F(3, 7) \times\\
& F(3, 3)\\


\end{aligned}

变成。如您所见，我们有很多共同点。因此，可以通过省略共同部分的计算来有效地计算它。
具体来说，$i = m - 1, m - 2, \dots ,0$ 是按顺序计算的，对于 $i = j$ 只计算 $F(U_{j + 1}, U_j)$ 并将结果相乘$i = j + 1$。

4.技巧3：避免多重精度操作

$F(L, U)$ 的计算是乘法，可能非常大。因此，通过设计乘法的顺序，尽可能避免繁重的多精度运算。

作为一个具体的例子，考虑 $F(11, 99)$ 的计算。
$F(11, 99)$ 大于 $64$ 位整数，因此将区间 $[11, 99)$ 分成两半，得到 $F(11, 55) \times F(55, 99)$。如果我们重复这个除法直到 $F(L, U)$ 适合 $64$ 位整数，

\begin{aligned}
F(11, 99) =& F(11, 33) \times F(33, 45) \times F(45, 55) \times F(55, 67) \times\\
&F(67, 77) \times F(77, 89) \times F(89, 99)
\end{aligned}

变成。
因此，在计算 $F(11, 99)$ 的 $43$ 乘法中，多精度运算是 $\boldsymbol{6}$次足够的。

从$11$依次相乘的方法中，$39$相乘时超过$64$位整数，所以$\boldsymbol{30}$次将对执行多精度算术运算。

现在我们知道了$odd\_part(n)$ 的计算算法。

5. two_exponent 的数学表示

$two\_exponent(n)$ 是 $n!$ 中 $2$ 的个数，所以这是 $(2 的倍数) + (4 的倍数) + \cdots$。

所以，

two\_exponent(n) = \sum_{i=1}^{m}{\left\lfloor\frac{n}{2^i} \right\rfloor}

是。 $m$ 是满足 $\left\lfloor\frac{n}{2^m}\right\rfloor \gt 0$ 的最大整数，所以 $\log{n}$ 就足够了。你可以计算 $two\ _exponent(n)$ 快。

但是，它可以更简洁地表达为

two\_exponent(n) = n - popcount(n)

是。但是，当 $n$ 以 $2$ 为基数表示时，$popcount(n)$ 是'1' 的位数。

$\boldsymbol{\sum_{i=1}^{m}{\left\lfloor\frac{n}{2^i} \right\rfloor} = n - popcount(n) 的证明}$

当 $n = 0$ 时很明显。
令 $n$ 为大于或等于 $1$ 的整数，$m$ 为最大整数，使得 $\left\lfloor\frac{n}{2^m}\right\rfloor \gt 0$。然后，使用长度为 $m + 1$ 的任意序列 $a = \{a_0, a_1, \cdots a_{m}\}$，由 $0, 1$ 组成，

n = \sum_{i=0}^{m}{a_{i} 2^{i}}

可以表示为（相当于以 2 美元为基础的符号。）

使用这个

\begin{aligned}
\sum_{i=1}^{m}{\left\lfloor\frac{n}{2^i} \right\rfloor} =& \sum_{i=1}^{m}{\left\lfloor\frac{ \sum_{j=0}^{m}{a_{j} 2^{j}}}{2^i} \right\rfloor}\\[3ex]
=& \sum_{i=1}^{m}{\left\lfloor\frac{ a_{0} 2^{0} + a_{1} 2^{1} + \cdots + a_{m} 2^{m}}{2^i} \right\rfloor}\\[1ex]\\
=& \sum_{i=1}^{m}\left\{{\left\lfloor\frac{ a_{0} 2^{0} + a_{1} 2^{1} + \cdots + a_{i - 1} 2^{i - 1}}{2^i} \right\rfloor} + \sum_{j = i}^{m}{a_{j}2^{j - i}}\right\}\\[3ex]
=& \sum_{i=1}^{m}{\sum_{j = i}^{m}{a_{j}2^{j - i}}}
\end{aligned}

变成。最后一行我们使用 $2^0 + 2^1 + \cdots + 2^{i-1} \lt 2^i$ 。

在这里，如果你考虑一下 $i 和 j$ 之和的范围，它看起来像下图中的红点。

Pythonで階乗 math.factorial のアルゴリズム

所以，

\begin{aligned}
\sum_{i=1}^{m}{\sum_{j = i}^{m}{a_{j}2^{j - i}}} =& \sum_{j=1}^{m}{\sum_{i = 1}^{j}{a_{j}2^{j - i}}}\\[1ex]
=& \sum_{j=1}^{m}{a_{j}2^{j} \sum_{i = 1}^{j}{2^{-i}}}
\end{aligned}

可以改写为
由于 $i$ 上的总和是几何级数的总和，

\begin{aligned}
\sum_{i = 1}^{j}{2^{-i}} =& \frac{\frac{1}{2}\left\{\left(\frac{1}{2}\right)^{j} - 1\right\}}{\frac{1}{2} - 1}\\[2ex]
=& 1 - 2^{-j}
\end{aligned}

是。所以

\begin{aligned}
\sum_{i=1}^{m}{\left\lfloor\frac{n}{2^i} \right\rfloor} =& \sum_{j=1}^{m}{a_{j}2^{j} (1 - 2^{-j})}
\end{aligned}

变成。此外，当 $j = 0$

a_{0}2^{0}(1 - 2^{-0}) = a_{0}\cdot1\cdot(1 - 1) = 0

所以你可以添加这个

\begin{aligned}
\sum_{i=1}^{m}{\left\lfloor\frac{n}{2^i} \right\rfloor} =& \sum_{j=0}^{m}{a_{j}2^{j} (1 - 2^{-j})}\\
=& \sum_{j=0}^{m}{a_{j}2^{j}} - \sum_{j=0}^{m}{a_{j}}
\end{aligned}

它可以得到。
$1$ 项是 $n$ 本身，$2$ 项是当 $n$ 以 $2$ 为基数表示时 '1' 的数量，所以这是 $popcount(n)$。

从以上

\sum_{i=1}^{m}{\left\lfloor\frac{n}{2^i} \right\rfloor} = n - popcount(n)

我说。

6.设备4：嵌入结果

如您所见，math.factorial() 的算法有点复杂。因此，当 $n$ 很小时，它比天真的计算要慢。另外，我认为当 $n$ 较小时，使用的机会更多。
作为针对这些的对策，如果 $n$ 很小，

SmallFactorials = [1, 1, 2, 6, 24, 120]

结果嵌入在源代码中，如例如，如果这调用math.factorial(4)

return SmallFactorials[4]

最好只是
math.factorial() 的实现嵌入了大约 $n \le 20$，阶乘结果适合 64 位整数。¹

七、总结

总结math.factorial(n)的算法计算$n!$的要点：

当 $n$ 小时使用嵌入结果（发明 4）
当 $n$ 很大时
- 除以奇数乘法和位移（设计 1）
- 由于同一区间的乘积多次出现，结果被重复使用（方案2）
- 使用分治法（设计 3）在尽可能避免多重精度运算的同时进行计算

8. 实施

在 Python 中实现它看起来像这样：²³

def F(L: int, U: int) -> int:
    """
    L, U: 奇数
    [L, U) の奇数の総積を分割統治法で計算する。
    """
    if L >= U: return 1

    max_bits = (U - 2).bit_length()  # 掛けられる最大の奇数のビット長
    num_operands = (U - L) // 2  # 掛けられる奇数の個数

    # [L, U) の奇数の総積のビット長は　max_bits * num_operands を超えない
    # これが long に収まれば多倍長演算を回避して計算できる
    if max_bits * num_operands < 63:  # 63 = sys.maxsize.bit_length()
        total = L
        for i in range(L + 2, U, 2):
            total *= i
        return total

    # 多倍長演算を回避するために分割して計算する
    mid = (L + num_operands) | 1
    left = F(L, mid)
    right = F(mid, U)

    return left * right

def calc_odd_part(n: int) -> int:
    """
    n! を (奇数) * (2の冪乗) と表したときの (奇数) の部分を計算する
    """

    result = 1
    L_i = 3
    tmp = 1  # F(3, U_i)
    m = n.bit_length() - 1  # n // (2 ** m) > 0 となる最大の整数

    for i in range(m - 1, -1, -1):
        # U_i は n//(2**i) より大きい最小の奇数
        U_i = ((n >> i) + 1) | 1

        # [1, U_i)　のうち、[1, L_i) は計算済みなので再利用し [L_i, U_i) のみ計算する
        tmp *= F(L_i, U_i)

        # 計算済みの範囲を更新 (L_{i} <- U_{i + 1})
        L_i = U_i

        result *= tmp

    return result


SmallFactorials = [1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800,
                   39916800, 479001600, 6227020800, 87178291200, 1307674368000,
                   20922789888000, 355687428096000, 6402373705728000,
                   121645100408832000, 2432902008176640000]


def factorial(n: int) -> int:
    """
    n! を計算する
    """
    # n が小さい場合は埋め込んだ結果を使う
    if n <= 20: return SmallFactorials[n]

    # n! = odd_part * (2 ** two_exponent) と表せる
    odd_part = calc_odd_part(n)
    two_exponent = n - bin(n).count('1')  # n - popcount(n)

    return odd_part << two_exponent



"""
test
"""
import math

for n in range(1000):
    assert factorial(n) == math.factorial(n)

9.速度对比

使用timeit 模块

factorial_naive()
factorial()
math.factorial()

比较速度使用的版本是 Python 3.8.2。

对于 $n = 10^1$ ，factorial()，math.factorial() 嵌入结果更快。
在 $n = 10^2$ 时，语言（Python 或 C）很重要，factorial_naive() 使用简单算法在同一个 Python 中更快。
当我们上升到 $n = 10^5$ 时，算法差异变大，语言差异与算法差异相比变得可以忽略不计。

Pythonで階乗 math.factorial のアルゴリズム

10. 结论

通过使用 Python 的math 模块检查计算阶乘的算法，我发现做了各种巧妙的设计。

如果解释，错误，建议等有任何错误，请告诉我。

正确的，它取决于执行环境中unsigned long类型的位数，但在本文中，unsigned long类型表示一个64位整数。↩
为了便于解释，它与math.factorial()的部分形状不同，但算法的关键点是相同的。另外，还有部分不符合Python中推荐的表示法，重点与目前的解释相匹配。↩
为了匹配math.factorial()，将F(L, U)的计算除以使得结果适合$2^{63}$，但是用Python编写时似乎有更好的除法。↩

原创声明：本文系作者授权爱码网发表，未经许可，不得转载;

原文地址：https://www.likecs.com/show-308627045.html