GIMPS梅森素数搜寻及相关算法综述

1. 梅森素数简介

梅森数是指形如 $2^n-1$ 的数，记为 $M_n$，如果一个梅森数是素数那么它被称为梅森素数（Mersenne prime）。

梅森数是根据17世纪法国数学家马林·梅森（Marin Mersenne）的名字命名的，他列出了 $n\le 257$ 的梅森素数，不过他错误地包括了不是梅森素数的$M_{67}$ 和 $M_{257}$，而遗漏了 $M_{61}$、$M_{89}$ 和 $M_{107}$。

现在我们更关注的是，梅森素数具有哪些重要性质。

当 $n$ 为合数时，$M_n$ 一定为合数（因为当 a 整除 b 时，$M_a$ 一定整除 $M_b$，反之亦然）。但当 $n$ 为素数时，$M_n$不一定皆为素数，比如 $M_2=2^2-1=3$ 和 $M_3=2^3-1=7$ 是素数，但 $M_{11}=2^{11}-1=2047=23\times 89$ 却不是素数。

另一方面，当 $M_n$ 为素数时，则 $n$ 一定为素数。基于此，为搜寻梅森素数，我们可以将目标限定在 $n$ 取素数的场合。

截至2022年8月，已知的梅森素数共有51个。其中最大的梅森素数是 $$2^{82589933}-1$$ 该数足有 24,862,048 位十进制数，也是已知最大的素数。从1997年至今，所有新的梅森素数都是由 互联网梅森素数大搜索(GIMPS) 分布式计算项目发现的。然而下面的结论也是十分有趣的：

设 $q\equiv 3mod4$ 为素数。则 $2q+1$ 是素数的充分必要条件是 $2q+1|M_q$ ，因此对于这些素数 $q$（除3以外），$M_q$ 不可能是素数。

前几个这样的素数 $q$ 为 11, 23, 83, 131, 179, 191, 239, 251, 359, 419, 431, 443, 491, 659, 683, 719, 743, 911, 1019, 1031, 1103, 1223, 1439, 1451, 1499, 1511, 1559, 1583, 1811, 1931, 2003, 2039, 2063, 2339, 2351, 2399, 2459, 2543, 2699, 2819, 2903, 2939, 2963, 3023, 3299, …

以上正是搜寻梅森素数时需要进一步排除的情形，同时也节省了很多计算时间，但这个长长的素数序列需要提前创建出来。关于梅森数的素性检验，目前最好的方法当属卢卡斯-莱默检验（Lucas–Lehmer primality test）。

该方法由爱德华·卢卡斯于1878年发现，并由德里克·亨利·莱默于1930年代作了改进，因此得名。该方法基于循环数列的计算，其原理是：

$M_n$ 为素数当且仅当 $M_n|S_{n-2}$ , 其中 $S_0=4，S_k=S_{k-1}^2-2，k>0$，此数列具体为 4, 14, 194, 37634, 1416317954, 2005956546822746114, …

例如，$M_3=7$，$S_1=14$，由于 $M_3|S_1$，因此 $M_3$ 为素数。看似简单的问题，随著素数 $p$ 值的增大，每一个梅森素数 $M_p$ 的产生都无比艰辛。然而，人们对于梅森素数搜寻仍乐此不疲。

2. GIMPS 项目

目前梅森素数搜寻最重要的途径当属 GIMPS 了，即因特网梅森素数大搜索。自1996年11月到现在，所有梅森素数皆由 GIMPS 发现。GIMPS 官网成立于1996年：

通过点击 Instructions 首先创建自己的账户：

账户创建好后，点击上图中的 Downloads 下载对应操作系统的 Prime95 软件，笔者使用的是2022年8月9日发布的 Prime95 version 30.8b16 版本。本软件无需安装，直接解压即可使用。

开启程序后如果没有自动分配计算任务，则点击 Advanced → Round off checking 即可开启第一个计算任务。

图中是当前作者的计算任务，正在通过 LL 算法完成梅森数 $M_{66384743}$ 的 Double-Check。至于Prime95 其它功能，可以查看 Help 帮助。

我们也可以加入一个计算团队，比如笔者加入了 GIMPSChina 计算团队。

如果搜寻到了新的梅森素数，那么你也可以领取奖励，有关奖励的最新说明如下：

当然，对于绝大多数爱好者，参与 GIMPS 搜寻梅森素数也绝不是为了领取奖金，因为发现一个新的梅森素数实在是太难了。截止2022年7月18日，所有低于 61427809 的测试都得到了验证，即 $M_n$ 中，$n\le 61,427,809$ 的梅森数都已经通过了再次验证，并于2021年10月6日确定出了第48个梅森素数 $M_{57885161}$。而后续三个数，目前还未确定是不是第49、50、51个梅森素数，因此在序号后面加上了红色星号，有关序号的确认，目前还在 GIMPS 的大量计算中。

从2018年至今将近四年过去了，仍没有发现 $M_{82589933}$ 之后的下一个梅森素数或者发现 $M_{82589933}$ 之前可能漏掉的梅森素数。截止2022年3月30日，所有低于1.09亿的测试至少进行了一次，即 $M_n$ 中，$n\le 109,460,867$ 的梅森数都已经测试了至少一次，仍未发现新的梅森素数。可见计算之艰巨，要发现一个新的梅森素数，当以年记。如果加入 GIMPS 搜寻项目的人越多，也越能加快新梅森素数的发现。

3. 卢卡斯-莱默测验

卢卡斯-莱默测验基于如下定理：

定理 设 $p$ 是一个奇素数，则梅森数 $M_p$ 是素数当且仅当 $M_p|S_{p-2}$，其中 $S_0=4，S_k=S_{k-1}^2-2，k>0$。

该定理的好处在于，将复杂的素性判断归结为简单的整除、迭代问题，因此非常高效地提升了素性检验的效率。如此好的方法，只对梅森素数有效，至于一般的大整数，其素性判断只能通过其他方法。比如概率方向的Miller-Rabin算法，2021年5月8日已测试出一个 8,177,207 位可能的素数 $R_{8177207}$。另一方面，椭圆曲线方向的ECPP算法，能证明一个大整数是素数，也是目前无条件证明大整数是素数的最有效的方法。2022年6月1日，一个 50001 位的素数 $10^{50000}+65859$ 已经得到了证明，也是目前通过 ECPP 证明的最大的素数。

以下是一段卢卡斯-莱默测验的伪代码：

Lucas_Lehmer_Test(p):
  s := 4;
  for i from 3 to p do s := s^2-2 mod 2^p-1;
  if s == 0 then
    2^p-1 is prime
  else
    2^p-1 is composite;

这个测试的理论是 Lucas 在19世纪70年代末提出的，然后在1930年左右由 Lehmer 提出了这个简单的测试。通过序列 $S_{p-2}$ 对 $M_p$ 取模，以节省大量运算时间。这个测试对于二进制计算机还非常理想，因为对 $M_p$ 的二进制除法可以只使用旋转和加法来完成。

下面我们给出卢卡斯-莱默测验的完整 Python 代码：

# 最后输出的 lucas 值等于 0 当且仅当 Mp 为素数
def lucas_lehmer(p):
    my_list=[4]
    value=2**p-1
    lucas=4
    for val in range(1, p - 1):
        lucas = ((lucas*lucas)-2) % value
        my_list.append(lucas)
    if lucas == 0:
        print("prime")
    else:
        print("composite")
    # print(my_list)
    # 仅打印最后的 lucas 值
    print(lucas) 

经测试，计算 lucas_lehmer(19937) 所花费的时间仅为 23.4672 秒，而这已是第24个梅森素数，该数首次发现于1971年3月4日，仅有 6,002 位十进制数。

为了进一步加快更大梅森素数的判定，我们将取模运算替换成位运算，完整代码如下：

# 更快速判定梅森素数
def lucas_lehmer_fast(n):
    m = 2**n - 1
    s = 4
    for i in range(2, n):
        sqr = s*s
        s = (sqr & m) + (sqr >> n)
        if s >= m:
            s -= m
        s -= 2
    return s == 0

经测试，同样计算 lucas_lehmer_fast(19937) 所花费的时间仅为 5.3317 秒，比之前的算法有了明显的改善。而计算 lucas_lehmer_fast(86243) 所花费的时间仅为 4分01.4661 秒，而这已是第28个梅森素数，该数首次发现于1982年9月25日，该素数足有 25,962 位！

以上，我们了解了梅森素数和 GIMPS 的基本知识，感受了自己编程验证梅森素数的过程，体会了随着 $M_p$ 中 $p$ 值的增大，相应的计算会越来越困难，个人的能力越来越有限，因此通过 GIMPS 加入梅森素数的搜寻工作将是一件十分不错的事情。

---

参考文献
[1]: List of Known Mersenne Prime Numbers
[2]: Mersenne Primes: History, Theorems and Lists