埃拉托斯散筛法算法程序

试编写一个程序，找出2->N之间的所有质数。希望用尽可能快的方法实现。

【问题分析】：

这个问题可以有两种解法：一种是用“筛子法”，另一种是“除余法”。

如果要了解“除余法”，请看另一篇文章《求质数之除余法(C语言描述)》。

这里我们来讨论一下用“筛法”来解决这个问题。

先来举个简单的例子来介绍一下“筛法”，求2~20的质数，它的做法是先把2~20

这些数一字排开：

2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

先取出数组中最小的数，是2，则判断2是质数，把后面2的倍数全部删掉。

2 |

3 5 7 9 11 13 15 17 19

接下来的最小数是3，取出，再删掉3的倍数

2 3 | 5 7 11 13 17 19

一直这样下去，直到结束。剩下的数都是素数。

筛法的原理是：

1.数字2是素数。

2.在数字K前，每找到一个素数，都会删除它的倍数，即以它为因子的整数。如果k 未被删除，就表示2->k-1都不是k的因子，那k自然就是素数了。

(1)除余法那篇文章里也介绍了，要找出一个数的因子，其实不需要检查2->k，只要从2->sqrt(k)，就可以了。所有，我们筛法里，其实只要筛到sqrt(n)就已经找出所有的素数了，其中n为要搜索的范围。

(2)另外，我们不难发现，每找到一个素数k，就一次删除2k, 3k, 4k,..., ik，不免还是有些浪费，因为2k已经在找到素数2的时候删除过了，3k已经在找到素数3的时候删除了。因此，当i＜k时，都已经被前面的素数删除过了，只有那些最小的质因子是k的那些数还未被删除过，所有，就可以直接从k*k开始删除。

(3)再有，所有的素数中，除了2以外，其他的都是奇数，那么，当i时奇数的时候，ik就是奇数，此时k*k+ik就是个偶数，偶数已经被2删除了，所有我们就可以以2k为单位删除步长，依次删除k*k, k*k+2k, k*k+4k, ...。

(4)我们都清楚，在前面一小段范围内，素数是比较集中的，比如1->100之间就有25个素数。越到后面就越稀疏。

因为这些素数本身值比较小，所以搜索范围内，大部分数都是它们的倍数，比如搜索1->100，这100个数。光是2的倍数就有50个，3的倍数有33个，5的倍数20个，7的倍数14个。我们只需搜索到7就可以，因此一共做删除操作50+33+20+14=117次，而2和3两个数就占了83次，这未免太浪费时间了。

所以我们考虑，能不能一开始就排除这些小素数的倍数，这里用2和3来做例子。

如果仅仅要排除2的倍数，数组里只保存奇数：1、3、5...，那数字k的坐标就是

k/2。

如果我们要同时排除2和3的倍数，因为2和3的最小公倍数是6，把数字按6来分

组：6n, 6n+1, 6n+2, 6n+3, 6n+4, 6n+5。其中6n, 6n+2, 6n+4是2的倍数，6n+3是3的倍数。所以数组里将只剩下6n+1和6n+5。n从0开始，数组里的数字就一次是1, 5, 7, 11, 13, 17...。

现在要解决的问题就是如何把数字k和它的坐标i对应起来。比如，给出数字89，它在数组中的下标是多少呢？不难发现，其实上面的序列，每两个为一组，具有相同的基数n，比如1和5，同是n=0那组数，6*0+1和6*0+5；31和35同是n=5那组，6*5+1和6*5+5。所以数字按6分组，每组2个数字，余数为5的数字在后，所以坐标需要加1。

所以89在第89/6=14组，坐标为14*2=28，又因为89%6==5，所以在所求的坐标上加1，即28+1=29，最终得到89的坐标i=29。同样，找到一个素数k后，也可以求出k*k的坐标等，就可以做筛法了。

这里，我们就需要用k做循环变量了，k从5开始，交替与2和4相加，即先是5+2=7，再是7+4=11，然后又是11+2=13...。这里我们可以再设一个变量gab，初始为4，每次做gab = 6 - gab，k += gab。让gab在2和4之间交替变化。另外，2和4都是2的幂，二进制分别为10和100，6的二进制位110，所以可以用k += gab ^= 6来代替。参考代码：

gab = 4;

for (k = 5; k * k <= N; k += gab ^= 6)

{

...

}

但我们一般都采用下标i从0->x的策略，如果用i而不用k，那应该怎么写呢？

由优化策略(1)可知，我们只要从k2开始筛选。n=i/2，我们知道了i对应的数字k

是素数后，根据(2)，那如何求得k2的坐标j呢？这里假设i为偶数，即k=6n+1。

k2 = (6n+1)*(6n+1) = 36n2 + 12n + 1，其中36n2+12n = 6(6n2+2n)是6的倍数，所以k2除6余1。

所以k2的坐标j = k2/6*2 = 12n2+4n。

由优化策略(2)可知，我们只要依次删除k2+2l×k, l = 0, 1, 2...。即

(6n+1)×(6n+1+2l)。

我们发现，但l=1, 4, 7...时，(6n+1+2l)是3的倍数，不在序列中。所以我们只要依次删除k2, k2+4l, k2+4l+2l...，又是依次替换2和4。

为了简便，我们可以一次就删除k2和k2+4l两项，然后步长增加6l。所以我们需要求len=4l和stp=6l。不过这里要注意一点，k2+4k=(6n+1)*(6n+5)，除以6的余数是5，坐标要加1。

len = k*(k+4)/6*2 - k2/6*2 = (6n+1)*(6n+1+4)/6*2+1 - (6n+1)*(6n+1)/6*2 = (12n2+12n+1) - (12n2+4n) = 8n+1;

stp = k*(k+6)/6*2 - k2/6*2 = 12n+2;