Python：为固定的 exp 和 mod 或矢量化加速 pow(base,exp,mod)

问题描述

我的代码的瓶颈是对非常大的整数重复调用 pow(base,exponent,modulus)（numpy 不支持这么大的整数，大约 100 到 256 位）。 但是，我的指数和模数总是相同的。我可以以某种方式利用它来加速自定义函数的计算吗？ 我试图定义一个如下所示的函数（下面的函数用于通用模数和指数）。

然而，即使我对每个操作进行硬编码，而没有针对固定指数和模数的 while 循环和 if 语句，它也比 pow 慢。

def modular_pow(self,base,modulus):
    result = 1
    base = base % modulus
    while exponent > 0:
        if (exponent % 2 == 1):
            result = (result * base) % modulus
        exponent = exponent >> 1
        base = (base * base) % modulus
    return result

另一种选择是，如果我能以某种方式“矢量化”它。我必须计算大约 100000000 个不同基值的 pow。虽然这些值经常在我的脚本运行之间发生变化（因此查找表将没有用），但我会在我运行的那一刻知道这些值（我可以一次计算它们）。

有什么想法吗？ 我通过使用 gmpy2 的 mpz 数据类型获得了一些加速，但它仍然太慢。

解决方法

好消息，坏消息。好消息是，当模数 m 固定时， 有多种方法可以加快计算 a*b % m。搜索“巴雷特还原”和“蒙哥马利还原”。它们以不同的方式工作，通过预先计算与 m 相关的常量，使得 % m 可以通过乘法和移位来计算，而无需除法。

坏消息：要找到余数，两种方法都需要（除了更便宜的操作）两次乘法。因此，除非乘法比除法便宜得多，否则他们不会支付整体费用。

出于这个原因，除非模数“真正”很大，否则它们通常会更慢-按照现代标准，您的“大约 100 到 256 位”仍然偏小，仅比原生 64 位机器宽几倍整数。诸如基于 FFT 的快速乘法之类的事情需要更大的整数才能得到回报。

CPython 的内置模块化 pow 已经在使用“二进制”方案，这与您在 Python 中编码的内容一致，但更高级（如果指数“足够大”，内置 pow 会将其视为在基数为 32，每次循环迭代消耗 5 个指数位）。

在 Python 中实现蒙哥马利约简并用蒙哥马利拼写替换代码中的模乘法的快速尝试中，modular_pow() 在模数增长到数万之前并没有比内置的快位。对于大约 256 位的输入，它大约慢了 3 倍。

这是一个混合包：Python 代码没有利用“base 32”技巧，这可以带来实质性的好处。但是对于足够大的输入，CPython 使用比朴素的 Karatsuba 乘法更快，无除法蒙哥马利拼写可以从中受益（CPython 的 int 除法无论输入大小都没有加速技巧，并且 CPython 的内置模块化 pow 总是使用除法来求余数）。

所以，简短的课程：我知道在 CPython 中没有任何明显的事情可以加速 pow(a,b,c) 的单个实例。可能有一些 C 编码的加密库有合适的东西，但我不知道。

但另一个好消息是您的问题是“令人尴尬的并行”。如果您有 N 个处理器，您可以为每个处理器提供 100000000/N 的输入，并且它们都可以全速并行运行。这将使速度提高大约 N 倍。

但坏消息是你的整数真的不是“大”（它们足够小，我敢打赌你仍然可以使用内置 pow 每秒计算数千个模块化 pow），并且进程间通信成本可能消除并行进行 N 次计算的好处。这完全取决于您如何获得输入以及您想对结果做什么。

跟进

《应用密码学手册》(HAC) chapter 14 基本上阐述了 gonzo 模幂算法的最新技术。

查看代码，GMP 已经实现了他们拥有的所有技巧。这包括我提到的事情（蒙哥马利约简，并使用大于 2 的 2 的幂来在每次循环迭代中咀嚼更多的指数位）。还有其他我没有提到的（例如，GMP 有一个特殊的内部平方程序，它可以节省可能不等整数的一般乘积的周期）。总之，这是一小堆实现代码。

我希望这就是您没有得到更多答案的原因：GMP 已经在做，在最坏的情况下，接近任何人想出的最佳方法。对您而言，加速并不是真正显着的，因为如前所述，您使用的整数实际上很小。

因此，如果您需要追求这一点，使用 GMP 可能是您获得的最快速度。如前所述，多处理是使用 N 个处理器获得理论 N 倍加速的一种显而易见的方法，但也注意到您没有说明上下文（这些输入来自何处或您需要对输出做什么）。因此，无法猜测这是否会为您带来回报。您需要的进程间通信越多，对潜在多处理加速的损害就越大。

注意：您所做的正是例如 RSA 公钥密码系统所做的，尽管它们通常使用更大的整数。也就是说，您的“基础”是他们的“消息”，而公共（或私有）RSA 密钥由固定指数和固定模数组成。只有基础（消息或加密位）因加密/解密实例而异。对于给定的键，指数和模数始终相同。

许多世界一流的数学家都研究过这个问题，世界一流的黑客为达到峰值速度编写了算法。这就是为什么你应该放弃希望 HAC 忘记提到的更快的方法;-)

投机

绘制到 RSA 的连接提醒我：RSA 解密在实践中不是以“显而易见”的方式进行。相反，私钥的持有者知道密钥模数的质因数分解（在 RSA 中，模数是两个不同但保密的大质数的乘积），并且可以用来显着加快对那个的取幂模数。

因此（无法猜测），如果您获得模数实例的方式可以有效地计算它们的质因数分解，那么当它们复合时，这可用于获得显着的加速。

不过，对于素数模数来说不是那么多。关于唯一高度潜在有价值的技巧是，对于p质数和a不能被p整除，

pow(a,p) == pow(a,b % (p-1),p)

如果 b 可以比 p 大得多，那可以节省无限的时间。之所以有效，是因为根据费马小定理，

pow(a,p-1,p) == 1

对于 p 质数和 a 不能被 p 整除。例如，

>>> p
2347
>>> assert all(p % i != 0 for i in range(2,p))  # that is,p is prime
>>> pow(345,1000000000000000000000000000000,p)
301
>>> 1000000000000000000000000000000 % (p-1)
1198
>>> pow(345,1198,p) # same thing,but much faster
301

对于复合模数，对其每个素功率因数进行大致相同的操作，然后通过中国剩余定理将结果粘贴在一起。

如果您认为您的问题可以利用这一点，请搜索“模幂中国余数”以找到一些好的说明。

,

看着 wiki page。您的实现似乎不正确。将这两个语句移出 else 显着提高了性能。

这是我从维基百科上得到的

def modular_pow(base,exponent,modulus):
    if modulus == 1:
        return 0
    else:
        result = 1
        base = base % modulus
        while exponent > 0:
            if exponent % 2 == 1:
                result = (result * base) % modulus
            exponent = exponent >> 1
            base = (base * base) % modulus
        return result

输出：

print(modular_pow(4,13,497))

445

,

可以使用加窗NAF方法预先计算a^2、a^3、...、a^(2^w-1)。现在，您有 n/w 轮产品，而不是 n 个产品和平方。例如，在 256 位的 modexp 中，当 w=4 时，我们进行 6 次预计算。但不是 256 个平方和乘积，我们有 256/4=64 个乘积。以 14 次预计算为代价，这是一笔可观的节省。现在 4 位是 2^4=16 个可能的值。但是 NAF 在 -w+1..w-1 范围内表示这些。注意负指数需要 a^(-1) 的模逆。因此，仅使用编码来缩放正值比额外的乘法或需要计算模逆更优化。注意a^0和a^1不需要计算。

某些预计算可能会根据 NAF 形式的指数进行优化，因为您事先知道需要哪些值。

w 值必须根据取幂而不是大小进行调整。最佳值可以根据 n/w vs. 2^w-1 计算，也可以凭经验确定。

我很惊讶问题的固定指数方面还没有得到解决。还有几篇关于这个确切主题的论文。椭圆曲线标量乘法中使用了相同的技术，尽管通常与基数相似的点是固定的，而不是此处与指数等效的标量。相同的方法适用于固定基数，但预计算可以离线完成并更有效地重复使用，而使用指数，它们是即时完成的。