|
发表于 2006-2-23 23:10:29
|
显示全部楼层
- 生成一个列表(Forming a list)
- --------------------------------------------------------------------------------
- 很容易证明,如果 2P-1 是素数,则 P 也一定是素数。因此,搜索梅森素数的第一步就是生成一个用于测试的素数指数列表。
- --------------------------------------------------------------------------------
- 试验分解因子(Trial Factoring)
- --------------------------------------------------------------------------------
- 下一步是通过寻找小因子来排除一些指数。有一个非常高效的算法判断一个数是否能整除 2P-1。例如,让我们看一下 47 是否能够整除 223-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
- 因此,223 = 1 mod 47。两边同时减 1,223-1 = 0 mod 47。因此我们知道 47 是一个因子,从而 223-1 不是素数。
- 可以证明梅森数有一个非常好的性质:2P-1 的任何因子 q 必定是 2kP+1 的形式,并且 q 除以 8 的余数一定是 1 或者 7。最后,一个高效的程序可以利用任何可能的因子 q 必须是素数这一事实。
- GIMPS 程序的分解因子代码使用修正的厄拉托森斯(Eratosthenes)筛法,利用一个二进位表示一个可能的 2kP+1 形式的因子。这个筛排除能够被大约 40,000 以下的素数整除的任何可能的因子。同样,表示除以 8 的余数是 3 或者 5 的可能的因子的二进位被清除。这个过程排除大约百分之九十五的可能的因子。剩下的可能的因子使用上面描述的高效的算法进行测试。
- 现在唯一的问题是要试验分解多少因子?答案取决于三个因素:分解因子的代价、发现一个因子的概率和素性测试的代价。我们使用以下公式:
- 分解因子的代价 < 发现因子的概率 * 2 * 素性测试的代价
- 也就是说,分解因子所花费的时间必须小于期望被节省的时间。如果能够发现一个因子,我们就能够避免进行首次素性测试和复查。
- 根据以前分解因子的数据,我们知道发现一个 2X 到 2X+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
- --------------------------------------------------------------------------------
- 用 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 = 3E*2*P。最后,检查 x-1 和 2P-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 的值。
- --------------------------------------------------------------------------------
- 卢卡斯-莱默测试(Lucas-Lehmer testing)
- --------------------------------------------------------------------------------
- 卢卡斯-莱默素性测试是非常简单的:如果 P > 2, 2P-1 是素数当且仅当 SP-2 = 0,其中,S0 = 4,SN = (SN-12 - 2) mod (2P-1)。例如,证明 27 - 1 是素数的过程如下:
- S0 = 4
- S1 = (4 * 4 - 2) mod 127 = 14
- S2 = (14 * 14 - 2) mod 127 = 67
- S3 = (67 * 67 - 2) mod 127 = 42
- S4 = (42 * 42 - 2) mod 127 = 111
- S5 = (111 * 111 - 2) mod 127 = 0
- 为了高效地实现卢卡斯-莱默测试,我们必须寻找对巨大的数进行平方及对 2P-1 取余的快速方法。自二十世纪六十年代后期以来,对巨大的数进行平方的最快速的算法是:把巨大的数分裂成小片形成一个大数组,然后执行快速傅里叶变换(FFT),逐项平方,然后再进行快速傅里叶逆变换(IFFT)。参见克努特的《计算机程序设计艺术》第二卷“乘法能有多快?”一节(译注:中文版第267页)。1994年1月,由理查德·克兰多尔(Richard Crandall)和巴里·费金(Barry Fagin) 合著的题为“离散加权变换和大整数算术”的计算数学文章,引入了无理底数 FFT 的概念。这个改进使得计算平方的速度提高两倍以上,允许使用较小的 FFT,并且这一过程中自动执行了对 2P-1 取余步骤。虽然由于英特尔公司的奔腾处理器体系结构的原因,GIMPS 程序使用浮点 FFT,但彼得·蒙哥马利(Peter Montgomery)给出的一个纯整数加权变换的方法也能够被使用。
- 正如上一段所提到的,GIMPS 使用汇编语言编写的浮点 FFT 算法,充分利用流水线和高速缓存。因为浮点运算是不精确的,在每次迭代后浮点值舍入到整数。本来该有的整数结果和程序计算出来的浮点结果之间的差异叫做“卷折误差”。如果卷折误差超过 0.5 则舍入将产生不正确的结果 - 这意味着必须使用更大的 FFT。GIMPS 程序的错误检查确保最大卷折误差不超过 0.4。不幸地,这种错误检查的代价相当高,以致于不能在每次平方后都进行检查。存在另外一种代价很低的错误检查。FFT 平方的一个性质是:
- (输入 FFT 值的和)2 = (输出 IFFT 值的和)
- 由于我们使用浮点数,我们必须将上式中的“等于”改为“约等于”。如果上式中两个值实质上不等,将给出一个在 readme.txt 文件中描述过的 SUMINP != SUMOUT 错误。如果输入 FFT 值的和是一个非法的浮点数(例如无穷大),将给出一个 ILLEGAL SUMOUT 错误。不幸地,这种错误检查无法发现我们将在下一节中描述的所有错误。
- 卢卡斯-莱默测试发现一个新的梅森素数的概率有多大?一个简单的估计是再次利用发现一个 2X 到 2X+1 之间的因子的概率大约是 1/X 的事实。例如,你已经使用试验分解因子证明 210000139-1 没有比 264 小的因子,那么它是素数的概率是: 没有 65 二进位因子的概率 * 没有 66 二进位因子的概率 * ... * 没有 5000070 二进位因子的概率,即:
- 64 65 5000069
- -- * -- * ... * -------
- 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”作为严重错误统计。
复制代码 |
|