問題描述

試編寫一個程式，找出 2→N 之間的所有質數（質數的概念請看這裡），用盡可能快的方法實現。

問題分析

這個問題可以有兩種解法：一種是用“篩子法”，另一種是從 2→N 逐一檢測出質數。如果要了解“除餘法”，請看另一篇文章《求質數 之 除餘法》。

先通過一個簡單的例子來介紹一下“篩法”，求 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 | 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)可知，我們只要從 k² 開始篩選。 n=i/2，我們知道了 i 對應的數字 k 是素數後，根據(2)，那如何求得 k² 的座標 j 呢？這裡假設 i 為偶數，即 k=6n+1。

k² = (6n+1)*(6n+1) = 36n² + 12n + 1，其中 36n²+12n = 6(6n²+2n) 是6的倍數，所以 k² 除 6 餘 1。

所以 k² 的座標 j = k²/6*2 = 12n²+4n。

由優化策略(2)可知，我們只要依次刪除 k²+2l×k, l = 0, 1, 2...。即 (6n+1)×(6n+1+2l)。

我們發現，但l=1, 4, 7...時，(6n+1+2l)是3的倍數，不在序列中。所以我們只要依次刪除 k², k²+4l, k²+4l+2l...，又是依次替換2和4。

為了簡便，我們可以一次就刪除 k² 和 k²+4l 兩項，然後步長增加6l。所以我們需要求 len=4l 和 stp=6l。不過這裡要注意一點，k²+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. 先列出所有的數。

02 03 04 05 06 07 08 09 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...。

02 03 04 05 06 07 08 09 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 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 17 19 23 25 29 31 35 37 41 43 47 49 53 55 59
02 03 05 | 07 11 13 17 19 23 29 31 37 41 43 47 49 53 59

5. 去掉5後，選出最小的數7，為素數，刪除7*7=49，...

02 03 05 | 07 11 13 17 19 23 29 31 37 41 43 47 49 53 59
02 03 05 | 07 11 13 17 19 23 29 31 37 41 43 47 53 59

6. 去掉7後，第一個數 11 的平方 121 大於 60，所以結束。剩下的數字全為素數。

02 03 05 07 11 13 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 <stdlib.h> #include <stdio.h> #include <time.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); return EXIT_SUCCESS; }

執行結果

5761455
Time: 0.300000

執行環境：Linux debian 2.6.26-1-686、GCC (Debian 4.3.2-1.1) 4.3.2、Intel Core 2、512MB RAM

演算法比較

現在，我們已經擁有初步改進的“篩法”和“除餘法”的函數了，把它們加到自己的函式庫裡。方便下次呼叫。這裡，我想說一下個人對這兩種演算法的使用經驗：

就時間效率上講，篩法絕對比除餘法高。比如上面的程式碼，可以在半秒內篩一億以內的所有素數。如果用除餘法來解決這樣的問題，絕對可以考驗一個人的耐性。因此，在搜尋空間比較大的時候，“篩法”無疑會是首選。

但篩法是以空間換時間，用除餘法，我們只要開一個可以容納結果的陣列就可以了，而篩法開的陣列要求可以容納整個搜尋範圍；另外，我們用“除餘法”得到的結果，是一個已經排好序的素數序列，如果要解決的問題需要用到這些連續的素數，而且搜尋範圍也不大，那顯然除餘法很適合。而“篩法”得到的結果，是一個布林型的表格，通過它，你可以很輕鬆的判斷某個數是不是素數，但如果你想知道這個素數的下一個素數是多大，可能要費點勁了。

聯絡方式

我的郵箱，歡迎來信（[email protected]）
我的Blogger（子清行）
我的Google Sites（子清行）
我的CSDN部落格（夢婷軒）
我的百度空間（夢婷軒）