线性筛法

 1 #include <iostream>  
 2 #include <cstring>  
 3 using namespace std;  
 4         
 5 bool isPrime[1000000];  
 6 int prime[1000000];  
 7 int MAX, total, cnt;  
 8         
 9 void makePrime() {  
10     memset(isPrime,true,sizeof(isPrime));  
11     memset(prime,0,sizeof(prime));  
12     total = cnt = 0;  
13     for(int i = 2; i <= MAX; ++i) {  
14         if(isPrime[i]) prime[total++] = i;  
15         for(int j = 0; j < total && i*prime[j] <= MAX; ++j) {  
16             ++cnt;  
17             isPrime[i*prime[j]] = false;  
18 //i此时不是素数,只是拓展用  
19             if(i%prime[j] == 0) break;  
20         }  
21     }  
22 }  
23         
24 int main(){  
25     while(cin>>MAX){  
26         makePrime();  
27         cout<<total<<' '<<cnt<<endl;  
28         //for(int i = 0; i < total; ++i) cout<<prime[i]<<endl;  
29     }  
30 }

View Code

上面为线性筛法的C++代码。运行可以发现输入输出满足total + cnt = MAXN - 1，说明外层循环执行MAXN次，内层循环cnt次，每个isPrime都只会被赋值一次false，完全O(n)的筛法。

整个算法最核心的一句是if(i%prime[j] == 0) break：当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已经不需要了。

 1 bool isPrime[MAXN];
 2 int label[MAXN], prime[MAXN];
 3 int n, total;
 4 
 5 void makePrime() {
 6     n = 100000;
 7     for(int i = 2; i <= n; ++i) {
 8         if(!label[i]) {
 9             prime[total++] = i;
10             label[i] = total;
11         }
12         for(int j = 0; j < label[i]; ++j) {
13             if(i * prime[j] > n) break;
14             label[i * prime[j]] = j + 1;
15         }
16     }
17     for(int i = 0; i < total; ++i) isPrime[prime[i]] = true;
18 }

View Code

新偷了一个线性筛法，多了一个数组，但是速度应该比上面那个快，因为不用求余了。原理一样，就是对任意数 x ，只被它最小的那个质因子剔除一次。

 1 void makePrime() {
 2     n = 10000000;
 3     int ed = sqrt(n + 0.5);
 4     memset(isPrime + 1, true, n * sizeof(bool));
 5     for(int i = 2; i <= n; ++i) if(isPrime[i]) {
 6         prime[total++] = i;
 7         if(i > ed) continue;
 8         for(int j = i * i; j <= n; j += i)
 9             isPrime[j] = false;
10     }
11 }

View Code

埃氏筛法，复杂度不明（据说是O(nloglogn)），比第二个那个稍微慢一点，不过跟第一个差不多（第一个常数大……）。

对于每一个素数p，只需要筛去p×p，p×(p+1)……就可以了，因为p×q(q＜p)已经在q的第一个素因子被找到的时候被筛去了。