第13章
埃拉托色尼筛法 – 计算素数表能有多快？
素数是除1与其自身外没有其它整除因子的自然数。素数在自然数集中的分布没有呈现出规律。正是这个不规则性许多世纪以来一直激发了数学家的好奇心与研究兴趣。
n以内的素数表列出自然数1到n范围内所有的素数，其开始部分如下：
2 3 5 7 11 13 17 19 23 29 31 37 41 ……
多年来人们提出了许多关于素数的问题，但是这些问题并未都得到了解决。这里有两个例子。
哥德巴赫(1694-1764)在1742年发表了他的一个有趣的发现:
“任意不小于4的偶整数可以分解为两个素数之和。”
例如，我们发现：
4=2+2， 6=3+3，8=3+5，10=3+7 等
这个命题要求对每个不小于4的偶整数至少给出一个相应的素数的和式，但实际上对于大多数偶数，这样的和式不止一个。下面的图是根据素数表画出来的，x-轴表示偶数，纵坐标反映了每个偶数对应的不同素数对的个数。
随着n增大，纵坐标值呈现略微递增的趋势。至今没有发现任何不小于4的偶数不满足此命题，但也没有能够证明这个命题对所有的不小于4的偶数一定成立。
高斯(1777-1855)通过对素数计数的方法研究素数的分布，并定义了如下函数：
p(n) = 1到n之间素数的个数
下图所示为此函数图像。
这样的函数称为阶梯函数，原因很明显。高斯定义了一个对任意大的n能尽可能逼近p(n)的连续函数。要描述高斯的方法并检查其结果是否有效，需要构造素数表。（这个问题已经解决了，但更深入的讨论超出了本书范围。）
今天，素数不仅仍然对数学家是个挑战，同时也有很大的应用价值。例如非常大的素数在电子密码学中发挥核心作用。
从思想到方法
就我们所知，一个名叫埃拉托色尼(公元前276-194)的古希腊人提出了第一个计算素数的算法。他是亚历山大城一位地位很高的学者，担任过历史上保存古代知识最为完备的亚历山大城图书馆的馆长。他与其他学者研究了那个时代天文、地理与数学方面最关键的问题，诸如地球的直径是多大？尼罗河发源于何处？怎样构造一个立方体使其体积为两个已知立方体体积之和，等等。
下面我们按照埃拉托色尼的做法将一个基本的思想转化为实用的方法。即使在他那个时代，人们也能在纸草或沙地上按这个方法进行计算。我们还要探究要计算一个足够大的素数表究竟能算得多快。所谓“大”是在109的数量级上，比方说10亿。
简单的思想
根据素数的定义，对于任意的非素数m，存在两个自然数i, k，满足：
2£i,k£m 且 i×k=m
根据这一事实，可以构造一个非常简单的素数表算法：
l 列出从2到n的所有整数；
l 算出从2到n的所有整数对i,k，算出其乘积；
l 将前面的整数列表中所有等于某个乘积的数删除。
显然这个算法就是我们要的素数算法，因为经删除能剩下的数不是任意i,k (2£i,k£n) 的乘积，因此它们一定是素数。
计算有多快？
为了对算法进行分析，我们更精确地写出算法的每一步，并按行排列如下：
素数表（基本版）
1 procedure PRIME NUMBER TABLE
2 begin
3 将从2到n的所有整数列在一个表中；
4 for i:=2 to n do
5 for k:=2 to n do
6 删除表中的数i×k;
7 endfor
8 endfor
9 end
在第6步中，如果数i×k在表中不存在，则不执行任何操作。
这个算法很容易在计算机上编程实现，我们便可以测量其执行时间。在linux PC(3.2GHz)运行得到如下的运行时间数据：
n 103 104 105 106
时间 0.00秒 0.20秒 19.4秒 1943.4秒
很明显，n增大10倍会导致时间开销大约增加100倍。这是可想而知的，因为当i和k各自得取值范围扩大10倍时，形成的乘积i×k的数量会增加100倍。
由此我们可以分析当n为109时所需的计算时间：我们必须在表中对应于106的时间值上再乘一个因子(109/106)2=106，结果是1943.4×106秒=61年7个月。显然在应用中这是不现实的。
算法将时间耗费在何处？
对{2,…,n} 中所有i,k计算i×k 仅对k³i计算i×k
图13.1 在一定范围内计算i×k
算法必须生成一定范围内所有的乘积。(见图13.1(a))
可实际上每个乘积i×k只需要一次。一旦从列表中删除，算法完全没有必要再生成同样的乘积。多余的工作出现在何处呢？当I,k的值恰好互换时就会出现多余的计算，例如对i=5, k=3计算了乘积，后面遇到i=3, k=5，因为乘法满足交换律，又计算出相等的乘积。因此，假如我们限定k³i就可以避免这样的冗余。(见图13.1(b))
这一想法立即可以省去一半代价，对计算素数表需要等上30年10个月仍然是太长了。在哪里还能省下更多时间呢？算法第6步什么情况下不会执行？当i×k>n时。列表中的数最大到n，所以不可能删除大于n的值。
因此算法的第5步中的k-循环只会在k满足i×k£n时才会执行，由此立即可知有效的k值范围是k£n/i。这也提示我们，同样可以限制i值的范围。同时考虑i£k£n/i可得i2£√n。当i足够大时，没有需要考虑的k值。需要考虑的i和k值的范围如下图所示：
现在算法可以改进为如下形式：
素数表(改进版)
1 procedure PRIME NUMBER TABLE
2 begin
3 将2到n的所有整数列在一个表中；
4 for i:=2 to ëÖnû do
5 for k:=2 to ën/iû do
6 从列表中删除数i×k
7 endfor
8 Endfor
9 End
( ë×û 表示下取整，因为I,k只能取整数值。)
算法快了多少呢? 改进的算法时间代价如下：
n 104 105 106 107 108 109
时间 0.00秒 0.01秒 0.01秒 2.3秒 32.7秒 452.9秒
收益还是客观的。我们可以仔细看看n=109的情况：只需等待7分半钟。让计算机去算吧，我们利用这段时间看看是否还能改进！
我们需要i取所有的值吗？
我们来更仔细地看一下算法第4行i-循环的执行情况:i保持一个固定值，而k遍历所有可能的取值(算法第5行)。在此过程中i×k依次为：
I2, i(i+1), i(i+2), i(i+3), ……
因此，当k-循环结束时，除i本身之外其它i的整数倍不可能还留在列表中。对于小于i的其它数的整数倍情况也一样，因为前面同样会被删除。
如果i不是素数。例如i=4:i×k的值包括16,20,24,…。由于4本身是2的整数倍，这些数也是2的整数倍。其实当i=4时这些数已经在前面被删除，不再需要做任何事。对于其他大于4的偶数情况也是如此。
再例如i=9:在i-循环中得到的乘积是9的整数倍，也必是3的整数倍。在i=3的循环中这些数已经被删除，这里再列出就多余了。这对于所有非素数都一样，他们所包含的质因子一定是前面已经遇到过的某个i值。由此可知，i-循环只需对素数执行。如下图所示：
至于i是否为素数可以在素数表中查，前提是表是完整的。那是否意味着我们首先必须完成表的计算才能信任它呢？
答案是“对也不对”。总地说是对的; 如果并不需要等算法结束那我们就可以简化算法。例如，当n=100,91不是素数，必须从列表中删除。只有当i=7,k=13时我们才能删除91，此时算法已进入最后一个i-循环了。
但我们现在面对是特殊情况。我们的目的并非判断任给的数是否为素数，只是判定特定的i是否为素数。而且我们需要判定的时刻也是特殊的:我们只需要在特定的i-循环开始执行k-循环之前判断i是否为素数。在此情况下，当前的素数表关于i是否为素数的信息一定是正确的。为什么？
从上面的讨论可知，对每个特定的i，删除的值满足i×k³i2。从另一个角度看，2,…,i2-1值域内的数不受影响。随着i值得增大，这个值域会扩大，并包含先前的值域。在下面的图中我们用蓝色标示这个值域。表的每行中第一个“错”数以红色标示。
所有蓝色标示的值域在算法结束前不再被改变。因此在当前i-循环开始执行k-循环之前，这些值域是正确的素数表。可以说在一个特定平方台阶上素数表完成了。我们要判定其素性的i数值用绿色标示。很显然它们一定在蓝色值域内。因此要判定i是否为素数只需要在当前的列表中查看。
由此我们可以改进算法中的i-循环：
素数表(埃拉托色尼版)
1 procedure PRIME NUMBER TABLE
2 begin
3 将2到n的所有整数列在一个表中；
4 for i:=2 to sqrt(n) do
5 if i在列表中 then
6 for k:=2 to ën/iû do
7 从列表中删除i×k
8 Endfor
9 Endif
10 Endfor
11 End
这就是那个聪明的希腊人给出的算法，称为“埃拉托色尼筛法”。叫“筛”是因为算法并不是直接“算出”素数，而是将非素数“筛掉”。
我们测量此算法的时间开销如下：
n 106 107 108 109
时间 0.02秒 0.43秒 5.4秒 66.5秒
即使n=109也只需要大约1分钟。
还能更快吗？
就像上面考虑如何限制i的值一样，我们也可以考虑如何限制k的取值:只需要考虑还在列表中的值。如果k是非素数，该从列表中删除，它必含质因子p<k。在i=p的i-循环中，p的所有整数倍，除其本省外，都已被删除。这其中一定包括k和k的倍数。因此没必要再考虑它们了。
很自然地我们会想到对算法的k-循环做如下的改变：
6 for k:=i to ën/iû do
7 if k在列表中 then
8 从列表中删除i×k
10 endif
但是，请当心！这想法有点误导。如果我们按照这个改变执行算法，会得到如下的结果：
2 3 5 7 8 11 12 13 17 19 20 23 27 28 29 31 32 37 ……
错在何处？让我们仔细查看一下算法执行开始的几步。首先我们建立了从2到n的整数列表(算法第3行)：
2 3 4 5 6 7 8 9 10 11……
第1步，设定i=2,接着k也被设置为2。2在列表中，因此我们会删除2×2=4:
2 3 - 4 5 6 7 8 9 10 11 ……
接下来k被设置为3。3也在列表中，因此我们会删除2×3=6：
2 3 - 5 - 7 8 9 10 11……
此时问题来了，4已被删除，不再列表中了。根据新的if的条件我们该跳过4去考虑k=5:
2 3 - 5 - 7 8 9 10 11……
这样8就被错误地保留下来，再也不会被删除了。问题在于一次k-循环只会删除i×k>k,然后k值递增。最后k的值会达到前面已经得到过的乘积i×k值(这个值在列表中已经不存在，因此它的倍数可能被漏删)。这方法给自身带来错误的影响。
解决方法是让k反向遍历其值域以避免它达到前面已经删除的值:
按照上述分析，可能考虑的乘积范围如下图所示：
综合起来可以得到以下的算法：
素数表(最终版)
1 procedure PRIME NUMBER TABLE
2 begin
3 将2到n的所有整数列在一个表中；
4 for i:=2 to sqlt(n) do
5 if i在列表中 then
6 for k:=sqlt(n/i) to i step -1 do
7 if k在列表中 then
8 从列表中删除i×k
9 endif
10 Endfor
11 Endif
12 Endfor
13 End
运行时间测试结果如下：
n 106 107 108 109
时间 0.01秒 0.15秒 1.6秒 17.6秒
这个结果按照今天的标准完全可以接受。以n=109为例，从最直观的方法开始，我们只是做了一些更深入的思考便将算法的速度提高了2.545亿倍。
我们能学到什么？
1. 简单的计算方法很可能效率不高。
2. 为了能提高速度，我们需要对采用的方法有更好的理解。
3. 往往存在几个不同的改进方法。
4. 数学思想会很有用处!
进一步的讨论
我们花费了一点时间对算法做了一点深入考虑，这17.6秒算得是个好结果。到底有多好呢？还能再进一步改进吗？
我们来梳理一下算法究竟做了什么。算法对不超过n的每个非素数至少生成一次并将其删除。在109以内的非素数准确地说有949,152,466个，而要计算的乘积的数量如下表所示：
基本版 改进版 埃拉托色尼版 最终版
乘积的数量 1018 9.44×109 2.55×109 9.49×108
与非素数数量的比 1.1×109 9.9 2.7 1.0
事实上，最终版执行的都是必须的计算，从这个意义上说这算法是最优的。
有趣的是这并非表明算法已不可能改进。非素数的数量并非绝对的度量标准，因为我们可以设法减少必须考虑的非素数数量。以下的技巧能发挥作用：开始的整数列表不必从2开始列出所有整数，只需包含2以及所有不小于3的奇数。既然我们知道不小于4的偶数不可能是素数，为什么还要将它们包括进来呢？在只包含奇素数候选对象的列表上执行算法，工作量
少很多。原来包含k-循环最多的i=2的循环完全可以省去了。现在只生成如下范围的乘积：
这个想法还可以进一步发挥，在开始的整数列表中也可以完全省略3的倍数(3本身除外)，开始的列表可以是如下的形式：
2 3 5 7 11 13 17 19 23 25 29 31 35 37 41 43 47 49 ……
“筛子”从i=5开始工作。
依照同样想法我们还可以在保持2到n的值域的前提下更进一步压缩初始列表的大小，比如省略5，7，11，13等等的倍数(自身除外)。
可以看出，随着列表的长度减小，计算时间与存储量也有相同的变化趋势:
省略
压缩 无压缩 2 2,3 2-5 2-7 2-11 2-13 2-17 2-19
运行时间(秒) 17.6 33.0 22.6 17.8 14.7 13.3 12.6 24.0 25.9
存储空间(兆字节) 119.2 59.6 39.7 31.8 27.3 24.8 22.9 21.6 20.4
这里用字位数组表示列表。数组中的位置i的值表示i的素性:表示素数，0则非素数。
如果将存储结构改为链表有两个缺点。一是要知道某个数的素性必须先找到该项；二是当n大到十进制9位时，我们要保存10亿以内的50,847,534个素数，占用空间高达1551.7兆字节。
还需注意一点，初始化时略去偶数会导致计算倍增。由于列表被压缩了，数组下标(1,2,3,4…)与实际考虑的整数(这里是2,3,5,7,11,…)之间必须转换。这要花些时间。不过总起来看，我们选择开始就省略2到13的倍数的做法还是产生了12.6秒这样非常好的结果。如果再进一步省略，转换的代价将超过表本身实际计算代价，那就没有意义了。
1万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
