Posted on 2010-08-07 17:03
MiYu 閱讀(1187)
評論(0) 編輯 收藏 引用 所屬分類:
ACM ( 數論 ) 、
ACM_資料
【問題描述】：
試編寫一個程序，找出2->N之間的所有質數。希望用盡可能快的方法實現。
【問題分析】：
這個問題可以有兩種解法：一種是用“篩子法”，另一種是“除余法”。
如果要了解“除余法”，請看另一篇文章《求質數 之 除余法(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的倍數全部刪掉。
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)可知，我們只要從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。
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.先列出所有的數。
2 3 4 5 6 7 8 9 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。
2 3 4 5 6 7 8 9 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 15 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 15 17 19 23 25 29 31 35 37 41 43 47 49 53 55 59
02 03 05 | 07 11 13 15 17 19 23 29 31 37 41 43 47 49 53 59
5.去掉5后，選出最小的數7，為素數，刪除7*7=49，
02 03 05 | 07 11 13 15 17 19 23 29 31 37 41 43 47 49 53 59
02 03 05 | 07 11 13 15 17 19 23 29 31 37 41 43 47 53 59
6.去掉7后，第一個數11的平方121大于60，所以結束。剩下的數字全為素數。
02 03 05 07 11 13 15 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 <time.h>
#include <stdio.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);
}
【運行結果】：
5761455
Time: 0.300000
運行環境：Linux debian 2.6.26-1-686、GCC (Debian 4.3.2-1.1) 4.3.2
【算法比較】：
現在，我們已經擁有初步改進的“篩法”和“除余法”的函數了，把它們加到自己的函數庫里。方便下次調用。
這里，我想說一下個人對這兩種算法的使用經驗：
就時間效率上講，篩法絕對比除余法高。比如上面的代碼，可以在半秒內篩一億以內的所有素數。如果用除余法來解決這樣的問題，絕對可以考驗一個人的耐性。因此，在搜索空間比較大的時候，“篩法”無疑會是首選。
但篩法是以空間換時間，用除余法，我們只要開一個可以容納結果的數組就可以了，而篩法開的數組要求可以容納整個搜索范圍；另外，我們用“除余法”得到的結果，是一個已經排好序的素數序列，如果要解決的問題需要用到這些連續的素數，而且搜索范圍也不大，那顯然除余法很適合。而“篩法”得到的結果，是一個布爾型的表格，通過它，你可以很輕松的判斷某個數是不是素數，但如果你想知道這個素數的下一個素數是多大，可能要費點勁了。
本文轉載自 :
版權聲明
本人的所有原創文章皆保留版權，請尊重原創作品。
轉載必須包含本聲明，保持本文完整，并以超鏈接形式注明原始作者“redraiment”和主站點上的本文原始地址。
聯系方式
我的郵箱，歡迎來信（redraiment@gmail.com）
我的Blogger（子清行）
我的Google Sites（子清行）
我的CSDN博客（夢婷軒）
我的百度空間（夢婷軒）