筛选法求质数

最近在做C语言名题精选一书中的习题，在书中题2.3、2.4、2.5中提到了有关质数的求解问题。在这里作出总结：

质数问题由来以久，有着许多解法，这里列出我所知道的解法。

 

1、传统的解法

 

传统的解法就是利用"一个数 n 如果是合数，那么它的所有的因子不超过sqrt(n)"这一性质来求解，时间复杂度为o(n*sqrt(n))，显然当n很大时这种方法便不太适用。

void getPrimeTrd(){
    int i,j;
    for(i = 2; i < MAX; ++i) {
　　　　　　 m = (int)sqrt(i);
        for(j = 2; j <= m; ++j){
            if(i % j == 0) {
                break;
　　　　　　　　  )
        }
        if(j > m) {
            printf("%d ",i);
           }
    }
}

 

习题2-3里则给出了上述算法的优化版。作出了以下优化

1、在外层循环除去了２和３的倍数，在任意６个数中只需要测试两个数据，如6,7,8,9,10,11其中6,8,10为２的倍数，9为3的倍数，只余下7和11。

2、在测试是否为质数时取prime[]中的数据，而除去2,3倍数的测试（外层循环就已经限制了不可能被2和3的倍数整除 ），只与prime[]中的质数相除。

void getPrimeTrdOpt(){
    int i,j,count;
    int prime[MAXSIZE];
    prime[0] = 2;
    prime[1] = 3;
    count = 2;
    for(i=5; i<MAXNUM; i+=2){
        if(i%3){
            for(j=0; prime[j]<=(int)sqrt(i); ++j){
                if(i%prime[j] == 0)
                    break;
            }
            if(prime[j] > (int)sqrt(i))
                prime[count++] = i;
        }
    }
    for(i=0; i<count; ++i)
        printf("%d ",prime[i]);
}

上面是本人根据题意的解法，下面是作者的解法。

 

void getPrimeTrdOpt2(){
     int  prime[MAXSIZE];   
     int  gap = 2;           
     int  count = 3;         
     int  may_be_prime = 5; 
     int  i;
     prime[0] = 2;           
     prime[1] = 3;          
     prime[2] = 5;
     while (prime[count] < MAXNUM) {
          may_be_prime += gap;
          gap = 6 - gap;
          for (i = 2; prime[i]*prime[i] <= may_be_prime; i++)
               if (may_be_prime % prime[i] == 0)
                    break;
          if (prime[i] > (int)sqrt(may_be_prime))     
               prime[count++] = may_be_prime; 
     }
     for (i = 0; i < count; i++) {
          printf("%d ", prime[i]);
     }   
}

显然作者的解答更加巧妙，其用了gap一个变量就达到了除去2和3倍数的功能（通过+2，+4，+2，+4....），观察除去2和3倍的连续数5,7,11,13,17.......,都是2,4间隔开的，当真十分巧妙。

2、筛选法求质数

筛选法又称筛法，是求不超过自然数N（N>1）的所有质数的一种方法。据说是古希腊的埃拉托斯特尼（Eratosthenes，约公元前274～194年）发明的，又称埃拉托斯特尼筛子。具体做法是：先把N个自然数按次序排列起来。1不是质数，也不是合数，要划去。第二个数2是质数留下来，而把2后面所有能被2整除的数都划去。2后面第一个没划去的数是3，把3留下，再把3后面所有能被3整除的数都划去。3后面第一个没划去的数是5，把5留下，再把5后面所有能被5整除的数都划去。这样一直做下去，就会把不超过N的全部合数都筛掉，留下的就是不超过N的全部质数。因为希腊人是把数写在涂腊的板上，每要划去一个数，就在上面记以小点，寻求质数的工作完毕后，这许多小点就像一个筛子，所以就把埃拉托斯特尼的方法叫做“埃拉托斯特尼筛”，简称“筛法”。（另一种解释是当时的数写在纸草上，每要划去一个数，就把这个数挖去，寻求质数的工作完毕后，这许多小洞就像一个筛子。）

 

版本一：

void getPrimeErato(){
    int prime[MAXSIZE];
    int i,j;
    for(i=0; i<MAXSIZE; ++i){
        if(i!=2 && i%2==0)
            prime[i] = NO;
        else
            prime[i] = YES;
    }
    for(i=3; i<=(int)sqrt(MAXSIZE); i+=2){
        if(prime[i]){
            for(j=i+i; j<MAXSIZE; j+=i){
                prime[j] = NO;
            }
        }
    }
    for(i=2; i<MAXSIZE; ++i){
        if(prime[i])
            printf("%d ",i);
    }
}

上面应该是最常见的筛选法求质数的程序，在初始化数组时将2的倍数设为了NO,这样做是为了简化程序。

 

版本二：

void getPrimeErato2(){
    int prime[MAXSIZE];
    int i,j;
    for(i=0; i<MAXSIZE; ++i){
        if(i!=2 && i%2==0)
            prime[i] = NO;
        else
            prime[i] = YES;
    }
    for(i=3; i<=(int)sqrt(MAXSIZE); i+=2){
        if(prime[i]){
            for(j=2; j<=MAXSIZE/i; ++j)
                prime[i*j] = NO;
        }
    }
 
    for(i=2; i<MAXSIZE; ++i){
        if(prime[i])
            printf("%d ",i);
    }
} 

版本二与版本一原理一样，版本二的可读性更好。

 

版本三（优化）：

void getPrimeEratoOpt(){
    int size = (MAXNUM - 3) / 2 + 1;
    int prime[size];
    int key;
    int i,j;
    for(i=0; i<size; ++i){
        prime[i] = YES;
    }
    for(i=0; i<size; ++i){
        if(prime[i]){
            key = i + i +3;
            for(j=key+2*i*i+4*i; j<size; j+=key){ 
                prime[j] = NO;
            }
        }
    }
    printf("%d ",2);   
    for(i=0; i<size; ++i){
        if(prime[i])
            printf("%d ",2*i+3);
    }
}

这里的程序，prime[i]不再表示i为质数式合数，而表示2i+3为质数式合数，如果prime[i] = OK，则2i+3为质数，否则为合数。这样做有什么好处？因为所在的质都是奇数，2i+3能表示所有的质数（2除外），这样做一方面减少了数组的空间，如需找出99以内的质数，用普通的筛选法则需定义一个100个长度的数组prime[100],而使用上面的程序只需定义一个（99-3）/2+1=49个长度的数组。空间个减少了一半。另一方面这做使得得到的质数全为奇数，直接排除了所有的偶数，提高了效率。程序中for(j=key+i;j<MAXSIZE/2-1;j+1=key){}这个for循环的作用是排除所有非质数，它是通过排除所有质数的奇数倍来排除奇数中的合数.具体操作如下，首先，确定一个质数2i+3,然后将2i+3的奇数倍对应的i值设为NO,如3（2i+3）为合数，则6i+9 = 2x+3,x=3i+3,于是将prme[3i+3]设为NO, 同理5（2i+3）=10+15=2x+3,x=5i+6，将prime[5i+6]设为NO。其中考虑到质数为5时会重复排除15（在3的时侯已经排除了），所以for循环中j=key+2*i*i+4*i,即从（2i+3）^2开始排除，而不是3（2i+3），为什么3 5 7 9 11 ... ...中的素数n的倍数（奇数）要从n开始？

就是说从1开始和从n开始的效果是一样，或者说n以前的奇数乘n都可以等于n以前（比n小的）素数的更大（比n大的）奇数倍数。由于n是奇数，可以设n以前（比n小的）的奇数为n-2i，比n大的奇数为n+2j；那么我们的任务就是，证明n*（n-2i）可以等于m*（n+2j），其中m为小于n的素数，i、j都是正整数。即n（n-m）=2mj+2ni。这有变成证

明n-m是偶数，这是显然的，两个奇数之差必然是偶数。

 

版本四（优化）：

void getPrimeEratoOpt2(){
    int i,j;
    int size = (MAXNUM-1)/2+1;
    int prime[size];
    for(i=0; i<size; ++i){
        prime[i] = YES;
    }
    for (i=3;i<size;i+=2)
    {
       if(prime[(i-1)/2]){
            for (j=i*i;j<MAXNUM;j+=i<<1)
                prime[(j-1)/2] = NO;
       }
    }
    printf("%d ", 2);
    for (i=1;i<size;++i)
    {
        if (prime[i])
            printf("%d ", (i<<1)+1);
    }
} 

版本四与版本三思路基本一致，区别在于一个用2i+3另一个用2i+1,最主要的区别在于在排除时过程与三正好相反，一个是prime[j] = NO,另一个是prime[(j-1)/2] = NO, 一个是循环直接控制i,而一个则是控制2i+1。从可读性上来说版本四的可读性更强，从空间利用率上来说版本三则更好。

 

版本五（线性筛选)

void getPrimeEratoLine(){
    int prime[MAXSIZE];
    int is_prime[MAXNUM];
    int count,i,j;
    count = 0;
    for(i=0; i<MAXNUM; ++i){
        is_prime[i] = 0;
    }
    for(i=2; i<MAXNUM ; ++i){
        if(is_prime[i] == 0){
            prime[count++] = i;
        }
        for(j=0; j<count && (i*prime[j]) < MAXNUM; ++j){
            is_prime[i*prime[j]] = 1;
            if(i % prime[j] == 0) break;
        }
    }
    for(i=0; i<count; ++i){
        printf("%d ",prime[i]);
    }   
}

这个算法有两层循环，第一层遍历2到n之间的所有自然数i，看看它是不是质数，如果是，则把i放进prime数组中。第二层循环是对所有未来的数进行筛选。对于当前正在处理的i，显然它乘以任何一个已经找到的素数的结果肯定是和数，它们将会被剔除。整个算法最核心的一句是第7行：当i可以被某个已经找到的质数整除的时候，循环退出，不再进行剔除工作。这样做的原因是：当prime[j]是i的因子的时候，设i=prime[j]*k，首先，我们可以肯定的说，prime[j]是i的最小质因数，这是因为第二层循环是从小到大遍历素数的；其次，我们可以肯定的说，i已经无需再去剔除prime[j']*i (j'>j) 形式的和数了，这是因为，prime[j']*i可以写成prime[j']*(prime[j]*k)=prime[j]*(prime[j']*k)，也就是说所有的prime[j']*i将会被将来的某个i'=prime[j']*k剔除掉，当前的i已经不需要了。
    然后我们再看看时间复杂度。虽然有两层循环，但是我们发现，正因为有了第7句，所有is_prime数组里面的所有false都只被赋值了一次，而且是在发现它的最小质因数的时候被赋值的。在外层循环执行了O(n)次操作的同时，内层循环里面总的操作次数也是O(n)次。因此，总的时间复杂度是O(n)。因为用了is_prime数组，空间复杂度也为O(n)。