问具有几百万位数的模数和指数的模幂的极快方法

EN

首先，re。答案1作者的评论“我不使用GMP，但我怀疑当他们使用FFT时，他们真的是指NTT”--不，当GMP说" FFT‘它意味着浮点FFT。IIRC也有一些基于NTT的例程，但对于bignum mul来说，它们与FFT没有竞争力。“

经过良好调整的FFT-mul比任何NTT都要高，因为由于舍入错误积累而导致的每字精度的轻微损失，远远超过现代CPU提供的高性能浮点能力的弥补，特别是当考虑到高性能的实现时，这些实现利用了x86_64族等CPU的向量-数学能力，而当前的迭代-- Intel Haswell、Broadwell和Skylake --具有大量的矢量浮点能力。(在这方面，我没有提到AMD，因为他们的AVX产品远远落后于英特尔( Intel )；他们的高潮期大约在2002年左右，从那以后，英特尔每年都在以越来越糟糕的方式击败他们。)GMP在这方面令人失望的原因是，相对来说，GMP的FFT是垃圾。我非常尊重GMP编码器的整体，但FFT的时间是FFT的时间，你不会得到分数的努力或例如有一个非常好的大加法。这里有一篇详细介绍GMP FFT-mul改进的论文：

Pierrick，Alex，Paul：“基于GMP的Sch nhage大整数乘法算法的实现”[http://www.loria.fr/~gaudry/publis/issac07.pdf]

这是从2007年开始的，但是AFAIK在下面的代码段中指出的性能差距并没有缩小；如果有的话，它已经扩大了。这篇论文很好地描述了可以使用的各种数学和算法改进，但让我们来谈谈金钱的报价：

“一个为整数乘法实现复杂浮点FFT的程序是George的Prime95。它主要是为了测试大型Mersenne数字2^p−1的原始性，在伟大的互联网默森素数搜索24。它使用by *2^n±c来表示乘法模a*2^n±c(a和c不太大)，参见图17。我们将Prime95版本24.14.2中的乘法模2^2wn−1与使用我们在奔腾4 ( 3.2 GHz )和Opteron 250 ( 2.4 GHz )上的SSA实现进行乘法的n字整数进行了比较，见图4。很明显，Prime95比我们的implementation要大得多，实际上在奔腾4上通常是10倍以上，在Opteron上则是2.5到3倍。“

接下来的几段是一堆节省面子的小说集。(我个人也认识3位作者中的2位，他们都是计算数论领域的佼佼者。)

请注意前面提到的George，他的Prime95代码从20年前问世不久就发现了所有的世界记录素数，他的核心的bignum例程已经以一个叫做GWNUM库的普通API化形式可用。您提到PFGW的速度比GMP快得多--这是因为PFGW使用GWNUM作为核心“重提升”算法，这就是PFGW中“GW”的来源。

我自己的FFT实现，它有泛型C构建支持，但就像乔治使用大量的x86向量-数学汇编程序的高性能的CPU家族，大约是60%-70%，比乔治的在目前的英特尔处理器系列。我相信这使得它成为x86上世界上第二快的bignum-mul代码。举个例子来说，我的代码目前正在使用一个30毫秒长的FFT (30*2^20倍)对一个大约2^29位的数字运行一个素数测试；因此，每个输入字超过17位。使用我的所有3.3 GHz哈斯韦尔4670四角的核心，它需要90毫秒/模。

顺便说一句，许多(如果不是大多数)世界顶尖的数学程序员都会在mersenneforum.org上闲逛，我鼓励你去看看，并向更广泛的(至少在这个特定领域)的专家们提问。我出现在与这里相同的手柄下；乔治·沃尔特曼( George Woltman )饰演“Prime95 95”，PFGW的马克·罗登基奇( Mark Rodenkirch )饰演“流氓”。

我根本不使用GMP，所以在处理这个问题时要记住这一点。

1. I更倾向于使用NTT代替FFT进行乘法。 它消除了舍入误差，并且与我的FFT实现相比，优化到相同点的实现更快。

- [C++ NTT implementation](https://stackoverflow.com/q/18577076/2521214)
- [C++ NTT sqr(mul) implementation](https://stackoverflow.com/q/18465326/2521214)

正如我提到的，我不使用GMP，但我怀疑当他们使用FFT时，他们实际上是指NTT (有限场傅立叶变换)。

1. modpow call.会导致您的测试和GMP原始测试的速度差异。 如果对它的调用太多，则会导致堆/堆栈崩溃，这会大大减慢速度。尤其是对bignums。尽量避免堆破坏，因此从操作数中删除尽可能多的数据，并对频繁调用的函数尽可能地返回调用。此外，有时还通过直接将函数源代码复制到代码中而不是调用(或者使用宏代替)来帮助消除瓶颈调用，只使用局部变量。 我认为GMP发布了它们的源代码，因此找到它们的modpow实现应该不会太难。你只需要正确地使用它
2. 只是为了澄清 您使用的是像20000+十进数字这样的大数字，这意味着每个数字的~8.4 KBytes。任何返回或非指针操作数都意味着将该数据量复制到堆堆栈或从堆堆栈中复制。这不仅需要时间，而且通常还会使CPU的缓存失效，这也会降低性能。 现在把它乘以算法迭代数，你就知道了。在调整许多bignum函数时也存在类似的问题，即使算法没有变化，加速比通常也比10000% (100次数)要高。只是限制/消除堆/堆栈垃圾。 因此，我认为您不需要更好的modpow实现，只需要更好地使用它，但是粗略地说，我可能错了，但是如果没有您正在使用的代码，就很难推断出更多。