[原創]求質數 之 篩法(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 | 3 5 7 9 11 13 15 17 19

接下來取出最小數 3，并删除 3 的倍數。

2 3 | 5 7 11 13 17 19

以此類推直至結束，剩餘之數皆為素數。

數字2是素數。

在數字k前，每找到一個素數，都會删除它的倍數，即以它為因子的整數。如果k未被删除，就表示2->k-1都不是k的因子，那k自然就是素數了。

除餘法那篇文章裡也介紹了，要找出一個數的因子，其實不需要檢查 2→k，隻要從 2->sqrt(k)，就可以了。所有，我們篩法裡，其實隻要篩到sqrt(n)就已經找出所有的素數了，其中n為要搜尋的範圍。

另外，我們不難發現，每找到一個素數 k，就一次删除 2k, 3k, 4k,..., ik，不免還是有些浪費，因為2k已經在找到素數2的時候删除過了，3k已經在找到素數3的時候删除了。是以，當 i＜k 時，都已經被前面的素數删除過了，隻有那些最小的質因子是k的那些數還未被删除過，所有，就可以直接從 k*k 開始删除。

再有，所有的素數中，除了 2 以外，其他的都是奇數，那麼，當 i 是奇數的時候，ik 就是奇數，此時 k*k+ik 就是個偶數，偶數已經被2删除了，所有我們就可以以2k為機關删除步長，依次删除 k*k, k*k+2k, k*k+4k, ...。

我們都清楚，在前面一小段範圍内，素數是比較集中的，比如 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。

最終，我們得到：

同理可以求出 k=6n+5 時的情況：

下面的代碼在實作上用了位運算，可能有點晦澀。

注：第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 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 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 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 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))

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;

運作環境：linux debian 2.6.26-1-686、gcc (debian 4.3.2-1.1) 4.3.2、intel core 2、512mb ram

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

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

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