求质数 之 筛法 ( 数论 C语言描述 zz )

【问题描述】：
　　　试编写一个程序，找出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 ;

　　　最终，我们得到：
　　　len  =  8n + 1 ;
　　　stp  =  12n + 2 ;
　　　 j　 =  12n2 + 4n;

　　　同理可以求出k = 6n + 5时的情况：
　　　len  =  4n + 3 ;
　　　stp  =  12n + 10 ;
　　　 j　 =  12n2 + 20n + 8 ;

　　　下面的代码在实现上用了位运算，可能有点晦涩。

★注：第5种优化方法还是理论阶段，下面的代码中并未采用这种优化算法，仅供大家参考。
　　　( 5 )由( 2 )可知，如果每找到一个素数k，能依次只删除以k为最小素数因子的数，那么每个数字就都只被删除一次，那这个筛法就能达到线性的O(n)效率了。比如数字600  =   2 * 2 * 3 * 5 * 11 ，其中2是它的最小素数因子。那这个数就被2删除了。 3 、 5 、11虽然都是它的因子，但不做删除它的操作。要实现这种策略，那每找到一个素数k，那从k开始，一次后面未被删除的数字来与k相乘，删除它们的积。比如要筛出2～60之间的素数：

　　　 1 .先列出所有的数。
　　　 2   3   4   5   6   7   8   9   10   11   12   13   14   15    50   51   52   53   54   55   56   57   58   59   60

　　　 2 .选出序列中的第一个数，即2，判断它是素数，然后从2开始，依次与剩下的未被删除的数相乘，删除它们的积。即2 * 2 = 4 ，  2 * 3 = 6 ， 2 * 4 = 8 。
　　　 2   3   4   5   6   7   8   9   10   11   12   13   14   15    50   51   52   53   54   55   56   57   58   59   60
　　　 02   |   03   05   07   09   11   13   15   17   19   21   23   25   27   29   31   33   35   37   39   41   43   45   47   49   51   53   55   57   59

　　　 3 .去掉2后，再选出序列中第一个数，即3，判断它是素数，然后从3开始，依次与剩下的数相乘，即3 * 3 = 9 ， 3 * 5 = 15 ， 3 * 7 = 21
　　　 02   |   03   05   07   09   11   13   15   17   19   21   23   25   27   29   31   33   35   37   39   41   43   45   47   49   51   53   55   57   59
　　　 02   03   |   05   07   11   13   15   17   19   23   25   29   31   35   37   41   43   47   49   53   55   59

　　　 4 .去掉3后，选出最小的数5，为素数，依次删除5 * 5 = 25 ， 5 * 7 = 35 ， 5 * 11 = 55 ，
　　　 02   03   |   05   07   11   13   15   17   19   23   25   29   31   35   37   41   43   47   49   53   55   59
　　　 02   03   05   |   07   11   13   15   17   19   23   29   31   37   41   43   47   49   53   59

　　　 5 .去掉5后，选出最小的数7，为素数，删除7 * 7 = 49 ，
　　　 02   03   05   |   07   11   13   15   17   19   23   29   31   37   41   43   47   49   53   59
　　　 02   03   05   |   07   11   13   15   17   19   23   29   31   37   41   43   47   53   59

　　　 6 .去掉7后，第一个数11的平方121大于60，所以结束。剩下的数字全为素数。
　　　 02   03   05   07   11   13   15   17   19   23   29   31   37   41   43   47   53   59   |

　　　上面的操作效率很高，但在计算机中模拟的时候却又很大的障碍：
　　　首先，计算机内存是一维的空间，很多时候我们不能随心所欲，要实现上面的算法，要求这个数据结构既能很高效地查找某个特定的值，又能不费太大代价对序列中的元素进行删除。高效地查找，用数组是最合适的了，能在O( 1 )的时间内对内存进行读写，但要删除序列中一个元素却要O(n)；单链表可以用O( 1 )的时间做删除操作，当然要查找就只能是O(n)了。所以这个数据结构很难找。
　　　其次，筛法的一个缺点就是空间浪费太大，典型的以空间换时间。如果我们对数组进行压缩，比如初始时就排除了所有偶数，数组0对应数字1，1对应3，。这样又会因为多了一道计算数字下标的工序而浪费时间。这又是一个矛盾的问题。
　　　也许我们可以试试折中的办法：数据结构综合数组和链表2种，数组用来做映射记录，链表来记录剩下的还未被删除的数据，而且开始也不必急着把链表里的节点释放掉，只要在数组里做个标记就可以了。下次遍历到这个数字时才删除。这样为了删除，可以算只遍历了一次链表，不过频繁地使用free()函数，也许又会减低效率。总之，我们所做的，依然是用空间来换时间，记录更多的信息，方便下次使用，减少再次生成信息所消耗的时间。

【程序清单】：

#include  < time.h >
#include  < stdio.h >

#define  N 100000000
#define  size (N/6*2 + (N%6 == 5? 2: (N%6>0)))

int  p[size  /   32   +   1 ]  =  { 1 };

int  creat_prime( void )
{
　　　　 int  i, j;
　　　　 int  len, stp;
　　　　 int  c  =  size  +   1 ;

　　　　 for  (i  =   1 ; ((i &~ 1 ) << 1 )  *  ((i &~ 1 )  +  (i >> 1 )  +   1 )  <  size; i ++ )
　　　　{
　　　　　　　　 if  (p[i  >>   5 ]  >>  (i  &   31 )  &   1 )  continue ;
　　　　　　　　len  =  (i  &   1 ) ?  ((i &~ 1 ) << 1 )  +   3 : ((i &~ 1 ) << 2 )  +   1 ;
　　　　　　　　stp  =  ((i &~ 1 ) << 1 )  +  ((i &~ 1 ) << 2 )  +  ((i  &   1 ) ?   10 :  2 );
　　　　　　　　j  =  ((i &~ 1 ) << 1 )  *  (((i &~ 1 ) >> 1 )  +  (i &~ 1 )  +   1 )  +  ((i  &   1 ) ?  ((i &~ 1 ) << 3 )  +   8   +  len: len);
　　　　　　　　 for  (; j  <  size; j  +=  stp)
　　　　　　　　{
　　　　　　　　　　　　 if  (p[j  >>   5 ]  >>  (j  &   31 )  &   1   ^   1 )
　　　　　　　　　　　　　　　　p[j  >>   5 ]  |=   1L   <<  (j  &   31 ),  -- c;
　　　　　　　　　　　　 if  (p[(j - len)  >>   5 ]  >>  ((j - len)  &   31 )  &   1   ^   1 )
　　　　　　　　　　　　　　　　p[(j - len)  >>   5 ]  |=   1L   <<  ((j - len)  &   31 ),  -- c;
　　　　　　　　}
　　　　　　　　 if  (j  -  len  <  size  &&  (p[(j - len)  >>   5 ]  >>  ((j - len)  &   31 )  &   1   ^   1 ))
　　　　　　　　　　　　p[(j - len)  >>   5 ]  |=   1L   <<  ((j - len)  &   31 ),  -- c;
　　　　}

　　　　 return  c;
}

int  main( void )
{
　　　　clock_t t  =  clock();

　　　　printf( " %d  " , creat_prime());
　　　　printf( " Time: %f  " ,  1.0   *  (clock()  -  t)  /  CLOCKS_PER_SEC);
}


【运行结果】：
5761455
Time:  0.300000

运行环境：Linux debian  2.6 . 26 - 1 - 686 、GCC (Debian  4.3 . 2 - 1.1 )  4.3 . 2

【算法比较】：
　　　现在，我们已经拥有初步改进的“筛法”和“除余法”的函数了，把它们加到自己的函数库里。方便下次调用。
　　　这里，我想说一下个人对这两种算法的使用经验：
　　　就时间效率上讲，筛法绝对比除余法高。比如上面的代码，可以在半秒内筛一亿以内的所有素数。如果用除余法来解决这样的问题，绝对可以考验一个人的耐性。因此，在搜索空间比较大的时候，“筛法”无疑会是首选。
　　　但筛法是以空间换时间，用除余法，我们只要开一个可以容纳结果的数组就可以了，而筛法开的数组要求可以容纳整个搜索范围；另外，我们用“除余法”得到的结果，是一个已经排好序的素数序列，如果要解决的问题需要用到这些连续的素数，而且搜索范围也不大，那显然除余法很适合。而“筛法”得到的结果，是一个布尔型的表格，通过它，你可以很轻松的判断某个数是不是素数，但如果你想知道这个素数的下一个素数是多大，可能要费点劲了。

本文转载自 :

---

版权声明

本人的所有原创文章皆保留版权，请尊重原创作品。
转载必须包含本声明，保持本文完整，并以超链接形式注明原始作者“redraiment”和主站点上的本文原始地址。

联系方式

我的邮箱，欢迎来信（redraiment@gmail.com）
我的Blogger（子清行）
我的Google Sites（子清行）
我的优快云博客（梦婷轩）
我的百度空间（梦婷轩）

转载于:https://www.cnblogs.com/MiYu/archive/2010/08/18/1802470.html