为什么 pow(a, d, n) 比 a**d % n 快那么多？

Question

我正在尝试实现Miller-Rabin素性测试，但对于中等大小的数字（约7位数），它花费的时间超过了20秒，让我感到困惑。最终，我发现以下代码行是问题的根源：

x = a**d % n

其中，a、d 和 n 都是相似但不相等的中等大小数，** 是指数运算符，% 是模运算符。

x = pow(a, d, n)

"

而相比之下，它几乎是瞬时的。

为了提供上下文，以下是原始函数：

"

from random import randint

def primalityTest(n, k):
    if n < 2:
        return False
    if n % 2 == 0:
        return False
    s = 0
    d = n - 1
    while d % 2 == 0:
        s += 1
        d >>= 1
    for i in range(k):
        rand = randint(2, n - 2)
        x = rand**d % n         # offending line
        if x == 1 or x == n - 1:
            continue
        for r in range(s):
            toReturn = True
            x = pow(x, 2, n)
            if x == 1:
                return False
            if x == n - 1:
                toReturn = False
                break
        if toReturn:
            return False
    return True

print(primalityTest(2700643,1))

一个计时的例子：

from timeit import timeit

a = 2505626
d = 1520321
n = 2700643

def testA():
    print(a**d % n)

def testB():
    print(pow(a, d, n))

print("time: %(time)fs" % {"time":timeit("testA()", setup="from __main__ import testA", number=1)})
print("time: %(time)fs" % {"time":timeit("testB()", setup="from __main__ import testB", number=1)})

输出结果（使用PyPy 1.9.0运行）：

2642565
time: 23.785543s
2642565
time: 0.000030s

输出（在Python 3.3.0下运行，2.7.2返回非常相似的时间）：

2642565
time: 14.426975s
2642565
time: 0.000021s

一个相关的问题是，为什么使用Python 2或3运行此计算几乎比使用PyPy快两倍，而通常情况下PyPy更快？

- lyallcooper

Answer 1

请看模数幂运算的维基百科文章。在进行 a**d % n 运算时，实际上需要计算非常大的数 a**d。但有一些方法可以计算a**d % n而无需先计算 a**d，这就是 pow 的作用。而 ** 运算符做不到这一点，因为它无法"看到未来"，即无法知道您会立即进行取模运算。

- BrenBarn

---

+1 这实际上就是文档字符串所暗示的意思：`>>> print pow.doc pow(x, y[, z]) -> number当只有两个参数时，等价于 xy。当有三个参数时，等价于 (xy) % z，但可能更高效（例如对于长整数）。` - Hedde van der Heide

根据你所使用的Python版本，这可能只在特定条件下成立。据我所知，在3.x和2.7中，你只能在整型（和非负幂）的情况下使用三个参数形式，并且你总是会得到具有本地“int”类型的模指数运算，但不一定适用于其他整型。但在早期版本中，存在适应C“long”的规则，允许使用三个参数形式对“float”等数据类型进行运算。 （希望你不要使用2.1或更早版本，并且不使用来自C模块的任何自定义整型类型，因此这些内容对你没有影响。） - abarnert

从你的回答中看来，似乎编译器无法看到表达式并进行优化，这并不正确。只是目前没有Python编译器这样做而已。 - danielkza

@danielkza：没错，我并不是想表明这在理论上是不可能的。也许用“不会预测未来”比“无法看到未来”更准确一些。请注意，优化在一般情况下可能非常困难甚至不可能。对于常数操作数，可以进行优化，但在x ** y % n中，x可能是一个实现了__pow__的对象，并且根据随机数返回多个实现了__mod__的对象，而这些对象以不同的方式取决于其他随机数等因素。 - BrenBarn

@danielkza：此外，这些函数的定义域不同：.3 ** .4％.5是完全合法的，但如果编译器将其转换为pow（.3，.4，.5），则会引发TypeError。编译器必须能够知道a，d和n保证是整数类型的值（或者可能只是特定类型的“int”，因为否则转换没有帮助），并且保证d为非负数。这是JIT可以想象做到的事情，但是对于动态类型且没有推断的语言的静态编译器来说，这是不可能的。 - abarnert

Answer 2

BrenBarn回答了您的主要问题。至于您的附言：

通常情况下，PyPy速度更快，但为什么在运行Python 2或3时却几乎快了一倍？

如果您阅读PyPy的性能页面，这正是PyPy不擅长的事情之一 - 实际上，他们给出的第一个示例就是：

糟糕的例子包括使用大型长整数进行计算-这是由无法优化的支持代码执行的。

从理论上讲，将一个巨大的指数幂和取模转换为模指数幂（至少在第一次操作后）可能是JIT可以进行的转换...但不是PyPy的JIT。

作为副产品，如果您需要使用巨大的整数进行计算，则可能希望查看第三方模块（如gmpy），它在某些非主流用途中比CPython的本地实现快得多，并且还具有许多其他功能，否则您可能需要自己编写，代价是不太方便。

- abarnert

---

长整型已经得到修复。尝试使用pypy 2.0 beta 1（它不会比CPython更快，但也不应该更慢）。gmpy没有处理MemoryError的方法 :( - fijal

@fijal：是的，gmpy 在某些情况下比原生 Python 更慢，并且使许多简单的事情变得不太方便。它并不总是最佳选择，但有时确实值得一试。因此，如果你正在处理大整数而 Python 的原生类型似乎不够快，那么它值得一看。 - abarnert

如果您不在意数字过大导致程序崩溃，那么请继续。 - fijal

@fijal：我不知道你为什么对gmpy这么愤怒。是的，如果我意外尝试计算Graham数，它最终会崩溃而不是引发可捕获但不可恢复的异常。我不敢想象在典型的64位CPU、8-32GB RAM和几TB交换空间的机器上，“最终”需要多长时间...但无论如何：有时这种差异很重要，有时则不重要。将其作为您决策的唯一因素，无论其是否实际相关，都是愚蠢的。 - abarnert

这是导致PyPy不使用GMP库的因素。这对你来说可能没问题，但对Python虚拟机开发人员来说不行。即使不使用大量RAM，malloc也可能失败，只需在那里放置一个非常大的数字。从那时起，GMP的行为是未定义的，而Python不能允许这种情况发生。 - fijal

@fijal：我完全同意不应该用它来实现Python内置类型。这并不意味着它永远都不应该被使用。 - abarnert

Answer 3

有一些用于进行模幂运算的快捷方式：例如，您可以为每个i从1到log(d)找到a**(2i) mod n，并将所需的中间结果（mod n）相乘。像具有3个参数的pow()这样的专用模幂函数可以利用这些技巧，因为它知道您正在进行模算术运算。Python解析器无法识别该裸表达式a**d % n，因此它将执行完整计算（这将需要更长时间）。

- atomicinf

Answer 4

计算 x = a ** d ％ n 的方式是将 a 的 d 次方计算出来，然后对 n 取模。首先，如果 a 很大，则会创建一个巨大的数字，然后被截断。然而，x = pow(a, d, n) 很可能被优化，只跟踪最后的 n 位数字，这些数字是计算取模运算所需的全部内容。

- Yuushi

---

计算 x 的 d 次方并不需要进行 d 次乘法运算，这种说法是错误的。实际上，只需要使用 O(log d)（非常少）次乘法运算即可完成计算。通过二分法进行幂运算可以在不使用模数的情况下实现。但是，大量的乘数才是导致计算时间增长的主要原因。 - John Dvorak

@JanDvorak 对，我不确定为什么我认为Python在**和pow中不会使用相同的指数算法。 - Yuushi

不是最后的“n”个数字，它只是在Z/nZ中进行计算。 - Thomas