Sieve of Eratosthenes

這是一個製作質數表的方法，聲名狼藉。

先將所有的正整數列出來。從2開始，將2的倍數全部移掉。接下來找下一個沒被移掉的數字，應該會找到3，將3的倍數全部移掉。接下來找下一個沒被移掉的數字，應該會找到5，將5的倍數全部移掉。接下來找下一個沒被移掉的數字，應該會找到7。將7的倍數全部移掉，如此不斷下去，最後剩下來的都是質數了。證明不難，不妨自己想想看。

bool prime[20000000]; void eratosthenes() { int i, j; for (i=0; i<20000000; i++) // 初始化 prime[i] = true; prime[0] = false; // 0 和 1 不是質數 prime[1] = false; for (i=2; i<20000000; i++) // 找下一個沒被移掉的數字 if (prime[i]) for (j=i+i; j<20000000; j+=i) // 將其倍數全部移掉 prime[j] = false; }

尚有個增進速度的方法。有個定理是這樣的：如果一個正整數x它不是質數，那麼x必定有個小於等於sqrt(x)的質因數。反過來說，只要將小於等於sqrt(x)的質因數找出來，並將這些質因數的倍數都移掉，則x一定是會被移掉的其中一個──也就是說，剩下來的數一定是質數。故篩選質數時，只需要找sqrt(20000000)以下的質數，將這些質數的倍數都刪掉，就可以得到20000000以內的所有質數了。

#include <cmath> // sqrt()函式在這裡面 bool prime[20000000]; void eratosthenes() { int i, j; for (i=0; i<20000000; i++) // 初始化 prime[i] = true; prime[0] = false; // 0 和 1 不是質數 prime[1] = false; int sqrt_value = sqrt(20000000); for (i=2; i<=sqrt_value; i++) // 找下一個沒被移掉的數字 if (prime[i]) for (j=i+i; j<20000000; j+=i) // 將其倍數全部移掉 prime[j] = false; }

時間複雜度是O(N*N)，空間複雜度是O(N)。

UVa 406 516 524 543 10140

Sieve of Eratosthenes

一般的篩法寫成這樣：

fill sieve[] with true; for (int i=2; i<MAX; i++) if (sieve[i]) for (int j=i+i; j<MAX; j+=i) sieve[j] = false;

接著我們一步步來改善並加速這個演算法。請比較前後程式碼的異同。

一、一個合數x必定只有小於等於sqrt(x)的質因數。換句話說，找出所有小於等於sqrt(MAX)的質因數，並篩掉他們所有的倍數之後，大於sqrt(MAX)且沒被篩掉的的數字必定為質數。

fill sieve[] with true; for (int i=2; i<sqrt_MAX; i++) if (sieve[i]) for (int j=i+i; j<MAX; j+=i) sieve[j] = false;

二、找到一個質數時，就要開始篩掉它所有的倍數，這是Sieve of Eratosthenes的精神。例如我們找到了11是質數，所以要刪掉11*2、11*3、11*4、……。但是，11的兩倍、三倍、四倍、……、十倍，早就應該在找到質數2、3、5、7的那時給篩掉了。也就是說，只要是比11小的倍率都早就篩掉了，所以直接從11的十一倍開始篩就可以了。

fill sieve[] with true; for (int i=2; i<sqrt_MAX; i++) if (sieve[i]) for (int j=i*i; j<MAX; j+=i) sieve[j] = false;

三、如果2的倍數都不做的話，那就可以省下處理2的倍數的時間。亦不需要再去砍質數的2n倍了。

fill sieve[] with true; for (int i=3; i<sqrt_MAX; i+=2) if (sieve[i]) for (int j=i*i; j<MAX; j+=i*2) sieve[j] = false;

四、既然2的倍數都不做，便可省下一半的記憶體空間了。令陣列的第0格代表數字1、第1格代表數字3、第2格代表數字5、……，以此類推。

// 這裡的 i, j 代表的是 index 而非數字 // 數字 = index * 2 + 1 index = (數字 - 1) / 2 fill sieve[] with true; for (int i=1; i < half_sqrt_MAX; i++) if (sieve[i]) for (int j=2*i*(i+1), k=2*i+1; j<half_MAX; j+=k) sieve[j] = false;

五、global變數會預設值為false。故乾脆將true和false倒過來用，這樣可以省去一開始全部要初始化為true的麻煩。

for (int i=1; i < half_sqrt_MAX; i++) if (!sieve[i]) for (int j=2*i*(i+1), k=2*i+1; j<half_MAX; j+=k) sieve[j] = true;

六、使用bit來取代bool陣列。一個unsigned int可以儲存32個bit，也就是說一個unsigned int可以當作32個欄位來使用。使用unsigned int陣列可將陣列的欄位縮小成原來的1/32，還可以節省記憶體空間和存取時間。

// sieve 是 unsigned int 陣列 // MAP(x) 可將原本 bool 陣列使用的 index 對應到正確的 bit 位置 #define MAP(x) ((sieve[x>>5]>>(x&31))&1) unsigned int i, j, k; for (i=1; i < half_sqrt_MAX; i++) if (MAP(i) == 0) for (j=2*i*(i+1), k=2*i+1; j<half_MAX; j+=k) sieve[j>>5] |= 1<<(j&31);

七、register關鍵字用來讓變數存放於CPU register，加速存取記憶體。

#define MAP(x) ((sieve[x>>5]>>(x&31))&1) register unsigned int i, j, k; for (i=1; i<half_sqrt_MAX; i++) if (MAP(i) == 0) for (j=2*i*(i+1), k=2*i+1; j<half_MAX; j+=k) sieve[j>>5] |= 1<<(j&31);

最後，要驗證一個數字是不是質數，可以這麼寫：

bool isprime(int p) { return p==2 || ( p>2 && (p&1)==1 && (MAP((p-1)>>1)==0) ) ; }

或者，在計算過程中，把算得的質數隨時儲存起來：

#define MAP(x) ((sieve[x>>5]>>(x&31))&1) register unsigned int i, j, k; vector<int> prime; for (i=1; i<half_sqrt_MAX; i++) if (MAP(i) == 0) { for (j=2*i*(i+1), k=2*i+1; j<half_MAX; j+=k) sieve[j>>5] |= 1<<(j&31); prime.push_back(i*2+1); } }

經過改進之後，時間複雜度變成O(NloglogN / logN)，空間複雜度仍為O(N)。

參考網址：http://www.algorithmist.com/index.php/Prime_Sieve_of_Eratosthenes

UVa 10311

6n±1 Method

這是一個精簡版的篩法，原理是：只拿2和3這兩個質數先篩過一遍，剩下的數字則用試除法驗證是不是質數。

2和3的最小公倍數是6，我們就把所有數字分為6n、6n+1、6n+2、6n+2、6n+3、6n+4、6n+5六種（n是倍率）。除以六會餘零的數字為6n，除以六會餘一的數字為6n+1，以此類推。

可以看出6n、6n+2、6n+3、6n+4會是2或3的倍數，不屬於質數。因此，只要驗證6n+1和6n+5是不是質數就可以了。（6n+5也可以寫成6n-1，意義不變。）

6n-1到6n+1，再到下一個6n-1，再到6n+1，把這些要驗證的數字由小排到大，可以發現之間的差值會是2 4 2 4 2 4 ...不斷跳二跳四。實作程式碼時，就可以直接用加法加二加四， 而不必用乘法及加減法計算6n-1、6n+1，如此一來程式的執行效率會好一點。

驗證的順序是：數字2和3明顯的是質數，不必驗證；若是從數字5開始驗證，那麼下一個要驗證的數字就是5+2，再下一個就是5+2+4，再下一個就是5+2+4+2，如此不斷下去。

// 以試除法判斷質數 bool isprime(int n) { for (int i=5; i*i<=n; i+=2) if (n%i == 0) return false; return true; } // 只檢查6n±1。這些數字的間隔為2 4 2 4.... void make_prime() { cout << "找到質數2"; cout << "找到質數3"; for (int i=5, gap=2; i<1000000; i+=gap, gap=6-gap) if (isprime(i)) cout << "找到質數" << i; } vector<int> prime; // 以試除法判斷質數 bool isprime(int n) { for (int i=0; prime[i]*prime[i]<=n; i++) if (n % prime[i] == 0) return false; return true; } // 只檢查6n±1。這些數字的間隔為2 4 2 4.... void make_prime() { cout << "找到質數2"; prime.push_back(2); cout << "找到質數3"; prime.push_back(3); for (int i=5, gap=2; i<1000000; i+=gap, gap=6-gap) if (isprime(i)) {cout << "找到質數" << i; prime.push_back(i);} }

這個方法的時間複雜度是O(NsqrtN)，空間複雜度是O(1)。事實上6n±1法比篩法慢上許多。不過6n±1法的程式碼構造較為單純，當要列舉的質數範圍不大時，有機會跑得比篩法還快。另外，6n±1法不需要開一條超大陣列來做計算，比起篩法它節省了很多空間，這也是它的優點。