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

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

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

                      删除最左 如果需要就     除以47
        平方          边二进位 乘以 2         的余数
        ------------  -------  -------------  ------
        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 是一个因子，从而 2²³-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        260
         5,160,000        261
         6,515,000        262
         8,250,000        263
        13,380,000        264
        17,850,000        265
        21,590,000        266
        28,130,000        267
        35,100,000        268
        44,150,000        269
        57,020,000        270
        71,000,000        271
        79,300,000        272

还有另外一个方法可被 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 程序使用该方法去寻找一些给人印象深刻的因子。例如：

        22944999-1 能够被 314584703073057080643101377 整除.
        314584703073057080643101377 等于 2 * 53409984701702289312 * 2944999 + 1.
        值 k, 53409984701702289312, 是非常平滑的:
        53409984701702289312 = 25 * 3 * 19 * 947 * 7187 * 62297 * 69061

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

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

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

评论