GIMPS:数学与算法

我们找到了新的梅森素数!

本页面讨论用于高效地搜索梅森素数的数学和计算机算法的一些知识。由于相对于数学家，我更多地是计算机程序员，因此我将不深入到太多的数学细节中，而是设法提供链接代替。

生成一个列表(Forming a list)

很容易证明，如果2^P-1是素数，则P也一定是素数。因此，搜索梅森素数的第一步就是生成一个用于测试的素数指数列表。

试验分解因子(Trial Factoring)

下一步是通过寻找小因子来排除一些指数。有一个非常高效的算法判断一个数是否能整除2^P-1。例如，让我们看一下47是否能够整除2²³-1。把指数23转换成二进制数，我们得到10111。从1开始，重复以下步骤：平方，删除指数的最左边二进位，如果该位是1，则将平方后得到的值乘以2，然后计算其除以47后的余数。

平方	删除最左边	二进位	如果需要就乘以 2	除以47的余数
1*1 = 1	1	0111	1*2 = 2	2
2*2 = 4	0	111	no	4
4*4 = 16	1	11	16*2 = 32	32
32*32 = 1024	1	1	1024*2 = 2048	27
27*27 = 729	1		729*2 = 1458	1

因此，2²³=1 mod 47。两边同时减1，2²³-1=0 mod 47。因此我们知道47是一个因子，从而223-1不是素数。

可以证明梅森数有一个非常好的性质：2^P-1 的任何因子_q必定是2kP+1的形式，并且_q除以8的余数一定是1或者7。最后，一个高效的程序可以利用任何可能的因子_q必须是素数这一事实。

GIMPS 程序的分解因子代码使用修正的厄拉托森斯(Eratosthenes)筛法，利用一个二进位表示一个可能的2kP+1形式的因子。这个筛排除能够被大约40,000以下的素数整除的任何可能的因子。同样，表示除以8的余数是3或者5的可能的因子的二进位被清除。这个过程排除大约百分之九十五的可能的因子。剩下的可能的因子使用上面描述的高效的算法进行测试。

现在唯一的问题是要试验分解多少因子？答案取决于三个因素：分解因子的代价、发现一个因子的概率和素性测试的代价。我们使用以下公式：

分解因子的代价<发现因子的概率*2*素性测试的代价

也就是说，分解因子所花费的时间必须小于期望被节省的时间。如果能够发现一个因子，我们就能够避免进行首次素性测试和复查。

根据以前分解因子的数据，我们知道发现一个2^X到2^X+1之间的因子的概率大约是1/X。本程序进行素性测试和分解因子所需的时间已经被计算出来。目前，本程序试图分解因子到：

指数上限	分解因子到
3,960,000	2⁶⁰
5,160,000	2⁶¹
6,515,000	2⁶²
8,250,000	2⁶³
13,380,000	2⁶⁴
17,850,000	2⁶⁵
21,590,000	2⁶⁶
28,130,000	2⁶⁷
35,100,000	2⁶⁸
44,150,000	2⁶⁹
57,020,000	2⁷⁰
71,000,000	2⁷¹
79,300,000	2⁷²

用P-1方法分解因子(P-1 Factoring)

还有另外一个方法可被GIMPS程序用来搜索因子，因而避免进行素性测试的花费。这个方法叫做波拉德(Pollard)(P-1)方法。如果_q是某数的一个因子，并且_q-1是高度复合的(也就是说_q-1只有小因子)，P-1 方法就可以找到因子_q。

该方法用于梅森数时甚至更高效。回忆一下，因子_q只能是2kP+1的形式。只要k是高度复合时，就很容易修改P-1方法去搜索因子_q。

P-1方法是十分简单的。在第一阶段我们挑选一个边界B1。只要k的所有因子都小于B1(我们称k为B1-平滑(B1-smooth))，P-1方法就能找到因子_q。我们首先计算E=(比B1小的所有素数的乘积)。然后计算x=3^E*2*P。最后，检查x-1和2^P-1的最大公约数，就可以知道是否找到一个因子。

使用第二个边界B2，我们可以改进波拉德算法，达到第二阶段。如果k在B1到B2之间刚好有一个因子，而其它因子都小于B1，我们就能够在第二阶段找到因子_q。这个阶段要使用大量的内存。

GIMPS程序使用该方法去寻找一些给人印象深刻的因子。例如：

2^2,944,999-1能够被314,584,703,073,057,080,643,101,377整除。

314,584,703,073,057,080,643,101,377等于2*53,409,984,701,702,289,312*2,944,999+1。

值k，53,409,984,701,702,289,312，是非常平滑的：

53,409,984,701,702,289,312=25*3*19*947*7,187*62,297*69,061

GIMPS如何智能地选择B1和B2呢？我们使用试验分解因子方法中的公式的变种。我们必须使下式取得最大值：

发现因子的概率*2*素性测试的代价-分解因子的代价

发现因子的概率和分解因子的代价都依赖于B1和B2的取值。当k是B1-平滑或者k是B1-平滑并且在B1到B2之间刚好有一个因子时，迪克曼(Dickman)函数(参见克努特(Knuth)的《计算机程序设计艺术》第二卷(译注:中文版第347页))用来确定发现因子的概率。本程序尝试许多B1的值，如果有足够的可用内存的话也尝试一些B2的值，用以确定使以上公式取得最大值的B1和B2的值。

卢卡斯-莱默测试(Lucas-Lehmer testing)

卢卡斯-莱默素性测试是非常简单的：如果P>2，2^P-1是素数当且仅当S_P-2=0，其中，S₀=4，S_N=(S_N-1²-2)mod(2^P-1)。例如，证明2⁷-1是素数的过程如下：

S₀=4

S₁=(4*4-2)mod 127=14

S₂=(14*14-2)mod 127=67

S₃=(67*67-2)mod 127=42

S₄=(42*42-2)mod 127=111

S₅=(111*111-2)mod 127=0

为了高效地实现卢卡斯-莱默测试，我们必须寻找对巨大的数进行平方及对2^P-1 取余的快速方法。自二十世纪六十年代后期以来，对巨大的数进行平方的最快速的算法是：把巨大的数分裂成小片形成一个大数组，然后执行快速傅里叶变换(FFT)，逐项平方，然后再进行快速傅里叶逆变换(IFFT)。参见克努特的《计算机程序设计艺术》第二卷“乘法能有多快？”一节(译注:中文版第267页)。1994年1月，由Crandall)和巴里·费金(Barry Fagin) 合著的题为“离散加权变换和大整数算术”的计算数学文章，引入了无理底数FFT的概念。这个改进使得计算平方的速度提高两倍以上，允许使用较小的FFT，并且这一过程中自动执行了对2^P-1取余步骤。虽然由于英特尔公司的奔腾处理器体系结构的原因，GIMPS程序使用浮点FFT，但彼得·蒙哥马利(Peter Montgomery)给出的一个纯整数加权变换的方法也能够被使用。

正如上一段所提到的，GIMPS使用汇编语言编写的浮点FFT算法，充分利用流水线和高速缓存。因为浮点运算是不精确的，在每次迭代后浮点值舍入到整数。本来该有的整数结果和程序计算出来的浮点结果之间的差异叫做“卷折误差”。如果卷折误差超过0.5则舍入将产生不正确的结果-这意味着必须使用更大的FFT。GIMPS程序的错误检查确保最大卷折误差不超过0.4。不幸地，这种错误检查的代价相当高，以致于不能在每次平方后都进行检查。存在另外一种代价很低的错误检查。FFT平方的一个性质是：

(输入FFT值的和)²=(输出IFFT值的和)

由于我们使用浮点数，我们必须将上式中的“等于”改为“约等于”。如果上式中两个值实质上不等，将给出一个在readme.txt文件中描述过的SUMINP!=SUMOUT错误。如果输入FFT值的和是一个非法的浮点数(例如无穷大)，将给出一个ILLEGAL SUMOUT错误。不幸地，这种错误检查无法发现我们将在下一节中描述的所有错误。

卢卡斯-莱默测试发现一个新的梅森素数的概率有多大？一个简单的估计是再次利用发现一个2^X到2^X+1之间的因子的概率大约是1/X的事实。例如，你已经使用试验分解因子证明2^10000139-1没有比2⁶⁴小的因子，那么它是素数的概率是：没有65二进位因子的概率*没有66二进位因子的概率*...*没有5000070二进位因子的概率，即：

化简后得到：64/5000070，或者1/78126。这个简单的估计不是很准确，它给出的公式是：(试验分解因子到多大的指数)/(指数/2)。进一步的工作表明更精确公式是：(试验分解因子到多大的指数-1)/(指数*欧拉常数(0.577...))。在上例中，是1/91623。这个更精确的公式是未经证明的。

复查(Double-checking)

为了核实首次的卢卡斯-莱默素性测试没有出错，GIMPS 程序运行第二次素性测试。在每次测试期间，最终的 SP-2 的最低 64 二进位，叫做余数，被打印出来。如果它们相同，GIMPS 宣称该指数已经被复查。如果它们不相同，素性测试被再次运行直到最后出现匹配。和首次测试相匹配的复查，通常是在首次测试之后大约两年进行。GIMPS 分配复查给较慢的计算机，因为该指数比正在进行的首次测试的指数小，以便较慢的计算机能够在合理的时间内完成其工作任务。

GIMPS 复查采取进一步的防护措施以避免程序设计错误。在开始卢卡斯-莱默测试之前，S0 的值被左移随机的二进位。每次平方刚好加倍我们左移的 S 值。注意对 2P-1 取余的步骤仅是简单地将第 P 位以上的位移到最低有效位，因此没有信息丢失。为什么我们要自找麻烦呢？因为如果计算 FFT 的程序代码有错误，对 S 值的随机的移位确保第二次素性测试中的 FFT 算法处理一个和首次素性测试完全不同的值。一个程序设计错误几乎不可能产生同样的最终 64 二进位余数。

历史上，卢卡斯-莱默测试运行期间没有报告严重错误时，结果的错误率大致是 1.5%。卢卡斯-莱默测试产生的错误被报告的比率大约是 50%。作为记录，我没有把“ILLEGAL SUMOUT”作为严重错误统计。