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() 如何计算阶乘。

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

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$ 之和的范围,它看起来像下图中的红点。

所以,

\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$ 时,算法差异变大,语言差异与算法差异相比变得可以忽略不计。

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

249人参与, 0条评论 登录后显示评论回复

你需要登录后才能评论 登录/ 注册