使用数百万位数的模和指数进行模幂运算的极快方法

Question

使用数百万位数的模和指数进行模幂运算的极快方法

performancealgorithmmathnumber-theorymodular-arithmetic

10

作为一个兴趣项目，我正在尝试寻找非常大的质数。这些质性测试包含模指数计算，即a^e mod n。为了简化说明，让我们称之为modpow操作。我想加速这个特定的计算。

目前我正在使用GMP的mpz_pown函数，但是它有点慢。我认为它太慢的原因是因为与软件PFGW进行完整的素性测试相比，对GMP的modpow的函数调用速度更慢（所以明确一下，这只是GMP的modpow部分，而不是我正在比较的整个自定义素性测试例程）。PFGW被认为是其领域中最快的，并且对于我的用例，它使用Brillhart-Lehmer-Selfridge素性测试-该测试也使用modpow过程-因此并非由于数学上的聪明才智使PFGW在这方面更快（如果我错了，请纠正我）。看起来GMP的瓶颈是modpow操作。对于具有略多于20,000位数字的数字，例如运行时间：GMP的modpow操作需要约45秒，而PFGW在9秒内完成了整个素性测试（涉及modpow）。随着数字变得更大，差异变得更加显著。GMP对于此测试比较使用FFT乘法和Montgomery减少，详见下文评论。

我做了一些研究。到目前为止，我了解到modpow算法使用平方求幂、整数乘法和模数规约——这些听起来都很熟悉。几个帮助方法可以提高整数乘法的运行时间：

为了改善平方求幂部分的运行时间，可以使用有符号数字表示法来减少乘法的次数（即将位表示为0、1或-1，并以这种方式表示位串，使其包含比其原始基数2表示中更多的零——这减少了平方求幂的运行时间）。

为了优化操作的模数部分，我知道以下方法：

蒙哥马利规约

所以这里是一个价值150,000美元的问题：是否有一个软件库能够高效地执行modpow运算，给定一个非常大的底数、指数和模数？（我希望能处理几百万位数）。如果您想提出建议，请尝试解释针对具有数百万位数的基数、模数和指数的情况的算法内部工作原理，因为一些库根据数字数量使用不同的算法。基本上，我正在寻找一个支持上述技术（或可能更聪明的技术）的库，它应该在运行算法时表现良好（至少比GMP好）。到目前为止，我已经搜索、发现并尝试了GMP和PFGW，但没有找到令人满意的解决方案（PFGW很快，但我只对modpow操作感兴趣，而且没有直接的编程接口）。我希望领域内的专家可以建议一个具备这些功能的库，因为似乎很少有库能够处理这些要求。

- webdevelopersdiary

mpn_powm 会使用一般的 mpn_mul 和 mpn_sqr，如果乘数的大小超过了 Toom-Cook 阈值，则会使用 FFT。除非 GMP 被编译时使用了 --disable-fft。mpn_powm 使用 Montgomery 模约简，包括偶模情况。因为 mpz_powm 会将模数分解为奇数因子和偶数 (2^i) 因子，然后对结果应用中国剩余定理。如果操作数很大，Montgomery 模约简例程 mpn_redc_n 也将使用 FFT。 - Brett Hale

感谢您的回复，Brett。看起来mpn_powm没有文档记录（如此提到）。该页面称mpn_powm比mpz_powm更快。您是否知道原因以及两者之间的区别？GMP函数索引也没有提到mpn_redc_n，听起来非常有用！我在哪里可以找到您提到的这些函数的文档？谢谢！ - webdevelopersdiary

我查了一下，这里提供参考：阈值在这里提到。在“最坏”的情况下，Toom-Cook的最大阈值似乎略小于600个limbs。一个limb可以容纳32或64位，因此这意味着它将在我的64位计算机上使用FFT来处理大于~ 600 * 64 = 39000位的整数。我还没有达到那么多位数。我将进行另一次测试以查看超过该限制会发生什么。 - webdevelopersdiary

我无法再编辑我的上一个评论，所以在新的评论中添加这个。我进行的约20,000位数字测试（在原帖中提到）是使用约70,000位数，因此我已经处于FFT频谱中了。 - webdevelopersdiary

mpz_ 函数是推荐的接口。mpn_ 函数执行底层例程，实现真正的工作。尽管如此，mpz_ 函数也设置工作内存、处理数据、为操作数选择最佳例程等。对于这样一个大型计算，绕过 MPZ 接口并不会节省任何东西，而且你还必须基本上复制设置代码。 - Brett Hale

我明白了，谢谢你的解释。现在我真的很好奇为什么GMP和PFGW的运行时间有这么大的差异！ - webdevelopersdiary

2个回答

4

请注意，我完全不使用GMP，因此请以此为前提考虑。

我更喜欢使用NTT而不是FFT进行乘法运算

它可以消除舍入误差，并且与我的FFT实现相比，在相同的优化点上可更快地完成计算。
- C++ NTT实现
- C++ NTT sqr(mul)实现
如我所述，我不使用GMP，但我怀疑当他们说他们使用FFT时，他们真正指的是NTT（有限域傅里叶变换）。
您测试的速度差异和GMP素数测试可能由于modpow调用引起。

如果对它的调用过多，则会导致堆栈崩溃，这会显着减慢速度。特别是对于bignums。尽量避免堆崩溃，因此对于经常调用的函数，请尽可能减少操作数和返回调用。有时，通过直接将函数源代码复制到您的代码中而不是调用（或使用宏）且仅使用局部变量来消除瓶颈调用可以帮助。

我认为GMP发布了他们的源代码，因此找到他们的modpow实现不应该太难。您只需要正确使用它
仅仅为了明确

您正在使用类似于20000+十进制数字的大数，这意味着每个数字的大小约为~8.4 KB。任何返回值或非指针操作数都需要将该数据从堆栈复制到内存中。这不仅耗费时间，而且通常会使CPU的CACHE无效，从而降低性能。

现在将此乘以算法迭代次数，您就会明白这一点。当我调整许多我的bignum函数时，也遇到过类似的问题，即使算法没有改变，速度提高的比例通常也超过10000％（100倍），只需限制/消除堆栈崩溃即可。

因此，我认为您不需要更好的modpow实现，而只需要更好地使用它，但是我可能错了，但是如果没有您正在使用的代码，则很难推断出更多信息。

- Spektre

@Webdevelopersdiary，是的，NTT是数字模板变换...与FFT相同，但作用于素数模数而不是复数域，这意味着数据量和操作量减少了2倍，但你需要n次单位根。在我的链接中也包含了我的C++实现类（在最开始优化得比我梦想的还要好）。我选择了最大的32位n次单位根，这消除了所有模除，并授予了最大的数据集使用权...此外，我强烈建议关注堆栈崩溃问题，如果你的代码浪费它，那么你可以获得更快的速度... - Spektre

你是否正在使用梅森素数#8，即2^31-1？我有一个问题：如果我们需要一个素数，这是否意味着我们只能使用世界上最大素数范围内的数字进行计算？（否则会出现溢出问题？）换句话说，它无法对比当前已知的（梅森？）素数更大的整数进行计算？ - webdevelopersdiary

1

我会谨慎地声称NTTs比FFT更快。因为实际上并非如此。即使是最好的NTTs，FFT也比它们快3-10倍，至少在尺寸达到内存成为问题之前。 - Mysticial

@Mysticial 已经进行了编辑，但是你有什么样的FFT优化想法（硬编码）？因为我的测量结果显示，我的NTT比我优化到同一点的FFT更快。所以要么FFT在专用硬件上运行，要么我错过了一些FFT优化...此外，FFT精度误差对于需要比NTT更精确的变量的大数非常大（我在12000+位上进行了测试，32位浮点数不够用），因此在我测试的每个方面中，NTT都获胜。我错过了什么？ - Spektre

哦，如果你正在使用32位浮点数，那么你完全走错了路。你需要使用双精度。在SSE2出现之前，使用80位扩展精度是标准做法。 - Mysticial

显示剩余11条评论

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- ewmayer · Accepted Answer

首先，关于答案1作者的评论“我不使用GMP，但我怀疑当他们说他们使用FFT时，他们真正指的是NTT” -- 不，当GMP说“FFT”时，它指的是浮点FFT。我IRC确实也有一些基于NTT的例程，但对于大整数乘法而言，这些例程与FFT相比效率较低。

经过精细调整的FFT-mul之所以能够击败任何NTT，原因在于由于舍入误差积累而导致的每个字精度的轻微损失被现代CPU提供的极其优越的浮点运算能力所弥补，特别是考虑到高性能实现所利用的诸如x86_64系列CPU的向量数学能力等方面，目前版本中的Intel Haswell、Broadwell和Skylake都具有极大的向量浮点能力。(我没有在这方面引用AMD，因为他们的AVX产品落后于Intel; 他们的高水平标志是2002年左右，自那时起，Intel每年都在逐渐恶化地超过他们。) GMP在这个领域表现失望的原因是GMP的FFT相对来说比较差。我非常尊重GMP编码人员的整体水平，但FFT计时是FFT计时，你不会因为努力或者拥有真正好的大数加法而得到分数。这里有一篇文章详细介绍了一系列GMP FFT-mul的改进：
Pierrick Gaudry, Alex Kruppa, Paul Zimmerman：“Schönhage-Strassen大整数乘法算法的基于GMP的实现” [http://www.loria.fr/~gaudry/publis/issac07.pdf]
这是2007年的内容，但据我所知，下面摘录的性能差距并没有缩小，如果有什么变化，那就是加剧了。该论文对可以部署的各种数学和算法改进进行了详细阐述，但让我们来看一下重要引用：一个实现复杂浮点FFT进行整数乘法的程序是George Woltman的Prime95。它主要是用于在Great Internet Mersenne Prime Search [24]中测试大型Mersenne数2^p - 1的素性。它使用DWT对模a*2^n ± c下的乘法进行计算，其中a和c不太大，详见[17]。我们比较了Prime95版本24.14.2中对2^2wn - 1取模的乘法和在Pentium 4 3.2 GHz和Opteron 250 2.4 GHz上使用我们的SSA实现的n字长整数的乘法，见图4。很明显，Prime95在Pentium 4上通常比我们的实现快一个数量级以上，在Opteron上则快2.5到3个数量级。

接下来的几段话是一连串的面子工程（还有一点，我个人认识这三位作者中的其中两位，他们都是计算数论领域的顶尖人物）。

请注意，前面提到的George Woltman，他的Prime95代码自推出以来已经发现了所有的世界纪录质数，并将其核心bignum例程以通用API方式提供，称为GWNUM库。你提到PFGW在FFT-mul方面比GMP快得多 - 这是因为PFGW使用GWNUM进行核心“重量级”算术运算，这就是PFGW中“GW”的含义所在。

我自己实现的FFT，具有通用-C构建支持，但像George的一样，在x86矢量数学汇编程序中使用大量的代码以获得高性能，目前在Intel处理器系列上比George的慢大约60-70％。我相信这使它成为x86上世界第二快的bignum-mul代码。例如，我的代码当前正在对具有大约2 ^ 29位的数字进行素性测试，使用30-Mdouble-length FFT（30 * 2 ^ 20个双精度浮点数）；因此每个输入字的位数略大于17位。使用我4个3.3 GHz Haswell 4670四核心的所有核心，每个modmul需要约90毫秒。

顺便说一下，许多（如果不是大多数）全球顶级的bignum-math编码人员都聚集在mersenneforum.org上，我鼓励您去看看并向那里更广泛的（至少在这个特定领域）专家观众提问。我在那里的用户名与此处相同；George Woltman出现为“Prime95”，PFGW的Mark Rodenkirch以“rogue”身份出现。