Posted on 2010-08-07 17:03
MiYu 閱讀(1178)
評(píng)論(0) 編輯 收藏 引用 所屬分類:
ACM ( 數(shù)論 ) 、
ACM_資料
【問題描述】：
試編寫一個(gè)程序，找出2->N之間的所有質(zhì)數(shù)。希望用盡可能快的方法實(shí)現(xiàn)。
【問題分析】：
這個(gè)問題可以有兩種解法：一種是用“篩子法”，另一種是“除余法”。
如果要了解“除余法”，請(qǐng)看另一篇文章《求質(zhì)數(shù) 之 除余法(C語言描述)》。
這里我們來討論一下用“篩法”來解決這個(gè)問題。
先來舉個(gè)簡(jiǎn)單的例子來介紹一下“篩法”，求2~20的質(zhì)數(shù)，它的做法是先把2~20這些數(shù)一字排開：
2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
先取出數(shù)組中最小的數(shù)，是2，則判斷2是質(zhì)數(shù)，把后面2的倍數(shù)全部刪掉。
2 | 3 5 7 9 11 13 15 17 19
接下來的最小數(shù)是3，取出，再刪掉3的倍數(shù)
2 3 | 5 7 11 13 17 19
一直這樣下去，直到結(jié)束。剩下的數(shù)都是素?cái)?shù)。
篩法的原理是：
1.數(shù)字2是素?cái)?shù)。
2.在數(shù)字K前，每找到一個(gè)素?cái)?shù)，都會(huì)刪除它的倍數(shù)，即以它為因子的整數(shù)。如果k未被刪除，就表示2->k-1都不是k的因子，那k自然就是素?cái)?shù)了。
(1)除余法那篇文章里也介紹了，要找出一個(gè)數(shù)的因子，其實(shí)不需要檢查2->k，只要從2->sqrt(k)，就可以了。所有，我們篩法里，其實(shí)只要篩到sqrt(n)就已經(jīng)找出所有的素?cái)?shù)了，其中n為要搜索的范圍。
(2)另外，我們不難發(fā)現(xiàn)，每找到一個(gè)素?cái)?shù)k，就一次刪除2k, 3k, 4k,
, ik，不免還是有些浪費(fèi)，因?yàn)?k已經(jīng)在找到素?cái)?shù)2的時(shí)候刪除過了，3k已經(jīng)在找到素?cái)?shù)3的時(shí)候刪除了。因此，當(dāng)i＜k時(shí)，都已經(jīng)被前面的素?cái)?shù)刪除過了，只有那些最小的質(zhì)因子是k的那些數(shù)還未被刪除過，所有，就可以直接從k*k開始刪除。
(3)再有，所有的素?cái)?shù)中，除了2以外，其他的都是奇數(shù)，那么，當(dāng)i時(shí)奇數(shù)的時(shí)候，ik就是奇數(shù)，此時(shí)k*k+ik就是個(gè)偶數(shù)，偶數(shù)已經(jīng)被2刪除了，所有我們就可以以2k為單位刪除步長，依次刪除k*k, k*k+2k, k*k+4k, 。
(4)我們都清楚，在前面一小段范圍內(nèi)，素?cái)?shù)是比較集中的，比如1->100之間就有25個(gè)素?cái)?shù)。越到后面就越稀疏。
因?yàn)檫@些素?cái)?shù)本身值比較小，所以搜索范圍內(nèi)，大部分?jǐn)?shù)都是它們的倍數(shù)，比如搜索1->100，這100個(gè)數(shù)。光是2的倍數(shù)就有50個(gè)，3的倍數(shù)有33個(gè)，5的倍數(shù)20個(gè)，7的倍數(shù)14個(gè)。我們只需搜索到7就可以，因此一共做刪除操作50+33+20+14=117次，而2和3兩個(gè)數(shù)就占了83次，這未免太浪費(fèi)時(shí)間了。
所以我們考慮，能不能一開始就排除這些小素?cái)?shù)的倍數(shù)，這里用2和3來做例子。
如果僅僅要排除2的倍數(shù)，數(shù)組里只保存奇數(shù)：1、3、5，那數(shù)字k的坐標(biāo)就是k/2。
如果我們要同時(shí)排除2和3的倍數(shù)，因?yàn)?和3的最小公倍數(shù)是6，把數(shù)字按6來分組：6n, 6n+1, 6n+2, 6n+3, 6n+4, 6n+5。其中6n, 6n+2, 6n+4是2的倍數(shù)，6n+3是3的倍數(shù)。所以數(shù)組里將只剩下6n+1和6n+5。n從0開始，數(shù)組里的數(shù)字就一次是1, 5, 7, 11, 13, 17。
現(xiàn)在要解決的問題就是如何把數(shù)字k和它的坐標(biāo)i對(duì)應(yīng)起來。比如，給出數(shù)字89，它在數(shù)組中的下標(biāo)是多少呢？不難發(fā)現(xiàn)，其實(shí)上面的序列，每兩個(gè)為一組，具有相同的基數(shù)n，比如1和5，同是n=0那組數(shù)，6*0+1和6*0+5；31和35同是n=5那組，6*5+1和6*5+5。所以數(shù)字按6分組，每組2個(gè)數(shù)字，余數(shù)為5的數(shù)字在后，所以坐標(biāo)需要加1。
所以89在第89/6=14組，坐標(biāo)為14*2=28，又因?yàn)?9%6==5，所以在所求的坐標(biāo)上加1，即28+1=29，最終得到89的坐標(biāo)i=29。同樣，找到一個(gè)素?cái)?shù)k后，也可以求出k*k的坐標(biāo)等，就可以做篩法了。
這里，我們就需要用k做循環(huán)變量了，k從5開始，交替與2和4相加，即先是5+2=7，再是7+4=11，然后又是11+2=13。這里我們可以再設(shè)一個(gè)變量gab，初始為4，每次做gab = 6 - gab，k += gab。讓gab在2和4之間交替變化。另外，2和4都是2的冪，二進(jìn)制分別為10和100，6的二進(jìn)制位110，所以可以用k += gab ^= 6來代替。參考代碼：
gab = 4;
for (k = 5; k * k <= N; k += gab ^= 6)
{
}
但我們一般都采用下標(biāo)i從0->x的策略，如果用i而不用k，那應(yīng)該怎么寫呢？
由優(yōu)化策略(1)可知，我們只要從k2開始篩選。n=i/2，我們知道了i對(duì)應(yīng)的數(shù)字k是素?cái)?shù)后，根據(jù)(2)，那如何求得k2的坐標(biāo)j呢？這里假設(shè)i為偶數(shù)，即k=6n+1。
k2 = (6n+1)*(6n+1) = 36n2 + 12n + 1，其中36n2+12n = 6(6n2+2n)是6的倍數(shù)，所以k2除6余1。
所以k2的坐標(biāo)j = k2/6*2 = 12n2+4n。
由優(yōu)化策略(2)可知，我們只要依次刪除k2+2l×k, l = 0, 1, 2。即(6n+1)×(6n+1+2l)。
我們發(fā)現(xiàn)，但l=1, 4, 7時(shí)，(6n+1+2l)是3的倍數(shù)，不在序列中。所以我們只要依次刪除k2, k2+4l, k2+4l+2l，又是依次替換2和4。
為了簡(jiǎn)便，我們可以一次就刪除k2和k2+4l兩項(xiàng)，然后步長增加6l。所以我們需要求len=4l和stp=6l。不過這里要注意一點(diǎn)，k2+4k=(6n+1)*(6n+5)，除以6的余數(shù)是5，坐標(biāo)要加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時(shí)的情況：
len = 4n+3;
stp = 12n+10;
j = 12n2+20n+8;
下面的代碼在實(shí)現(xiàn)上用了位運(yùn)算，可能有點(diǎn)晦澀。
★注：第5種優(yōu)化方法還是理論階段，下面的代碼中并未采用這種優(yōu)化算法，僅供大家參考。
(5)由(2)可知，如果每找到一個(gè)素?cái)?shù)k，能依次只刪除以k為最小素?cái)?shù)因子的數(shù)，那么每個(gè)數(shù)字就都只被刪除一次，那這個(gè)篩法就能達(dá)到線性的O(n)效率了。比如數(shù)字600 = 2*2*3*5*11，其中2是它的最小素?cái)?shù)因子。那這個(gè)數(shù)就被2刪除了。3、5、11雖然都是它的因子，但不做刪除它的操作。要實(shí)現(xiàn)這種策略，那每找到一個(gè)素?cái)?shù)k，那從k開始，一次后面未被刪除的數(shù)字來與k相乘，刪除它們的積。比如要篩出2～60之間的素?cái)?shù)：
1.先列出所有的數(shù)。
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.選出序列中的第一個(gè)數(shù)，即2，判斷它是素?cái)?shù)，然后從2開始，依次與剩下的未被刪除的數(shù)相乘，刪除它們的積。即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后，再選出序列中第一個(gè)數(shù)，即3，判斷它是素?cái)?shù)，然后從3開始，依次與剩下的數(shù)相乘，即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后，選出最小的數(shù)5，為素?cái)?shù)，依次刪除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后，選出最小的數(shù)7，為素?cái)?shù)，刪除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后，第一個(gè)數(shù)11的平方121大于60，所以結(jié)束。剩下的數(shù)字全為素?cái)?shù)。
02 03 05 07 11 13 15 17 19 23 29 31 37 41 43 47 53 59 |
上面的操作效率很高，但在計(jì)算機(jī)中模擬的時(shí)候卻又很大的障礙：
首先，計(jì)算機(jī)內(nèi)存是一維的空間，很多時(shí)候我們不能隨心所欲，要實(shí)現(xiàn)上面的算法，要求這個(gè)數(shù)據(jù)結(jié)構(gòu)既能很高效地查找某個(gè)特定的值，又能不費(fèi)太大代價(jià)對(duì)序列中的元素進(jìn)行刪除。高效地查找，用數(shù)組是最合適的了，能在O(1)的時(shí)間內(nèi)對(duì)內(nèi)存進(jìn)行讀寫，但要?jiǎng)h除序列中一個(gè)元素卻要O(n)；單鏈表可以用O(1)的時(shí)間做刪除操作，當(dāng)然要查找就只能是O(n)了。所以這個(gè)數(shù)據(jù)結(jié)構(gòu)很難找。
其次，篩法的一個(gè)缺點(diǎn)就是空間浪費(fèi)太大，典型的以空間換時(shí)間。如果我們對(duì)數(shù)組進(jìn)行壓縮，比如初始時(shí)就排除了所有偶數(shù)，數(shù)組0對(duì)應(yīng)數(shù)字1，1對(duì)應(yīng)3，。這樣又會(huì)因?yàn)槎嗔艘坏烙?jì)算數(shù)字下標(biāo)的工序而浪費(fèi)時(shí)間。這又是一個(gè)矛盾的問題。
也許我們可以試試折中的辦法：數(shù)據(jù)結(jié)構(gòu)綜合數(shù)組和鏈表2種，數(shù)組用來做映射記錄，鏈表來記錄剩下的還未被刪除的數(shù)據(jù)，而且開始也不必急著把鏈表里的節(jié)點(diǎn)釋放掉，只要在數(shù)組里做個(gè)標(biāo)記就可以了。下次遍歷到這個(gè)數(shù)字時(shí)才刪除。這樣為了刪除，可以算只遍歷了一次鏈表，不過頻繁地使用free()函數(shù)，也許又會(huì)減低效率。總之，我們所做的，依然是用空間來換時(shí)間，記錄更多的信息，方便下次使用，減少再次生成信息所消耗的時(shí)間。
【程序清單】：
#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);
}
【運(yùn)行結(jié)果】：
5761455
Time: 0.300000
運(yùn)行環(huán)境：Linux debian 2.6.26-1-686、GCC (Debian 4.3.2-1.1) 4.3.2
【算法比較】：
現(xiàn)在，我們已經(jīng)擁有初步改進(jìn)的“篩法”和“除余法”的函數(shù)了，把它們加到自己的函數(shù)庫里。方便下次調(diào)用。
這里，我想說一下個(gè)人對(duì)這兩種算法的使用經(jīng)驗(yàn)：
就時(shí)間效率上講，篩法絕對(duì)比除余法高。比如上面的代碼，可以在半秒內(nèi)篩一億以內(nèi)的所有素?cái)?shù)。如果用除余法來解決這樣的問題，絕對(duì)可以考驗(yàn)一個(gè)人的耐性。因此，在搜索空間比較大的時(shí)候，“篩法”無疑會(huì)是首選。
但篩法是以空間換時(shí)間，用除余法，我們只要開一個(gè)可以容納結(jié)果的數(shù)組就可以了，而篩法開的數(shù)組要求可以容納整個(gè)搜索范圍；另外，我們用“除余法”得到的結(jié)果，是一個(gè)已經(jīng)排好序的素?cái)?shù)序列，如果要解決的問題需要用到這些連續(xù)的素?cái)?shù)，而且搜索范圍也不大，那顯然除余法很適合。而“篩法”得到的結(jié)果，是一個(gè)布爾型的表格，通過它，你可以很輕松的判斷某個(gè)數(shù)是不是素?cái)?shù)，但如果你想知道這個(gè)素?cái)?shù)的下一個(gè)素?cái)?shù)是多大，可能要費(fèi)點(diǎn)勁了。
本文轉(zhuǎn)載自 :
版權(quán)聲明
本人的所有原創(chuàng)文章皆保留版權(quán)，請(qǐng)尊重原創(chuàng)作品。
轉(zhuǎn)載必須包含本聲明，保持本文完整，并以超鏈接形式注明原始作者“redraiment”和主站點(diǎn)上的本文原始地址。
聯(lián)系方式
我的郵箱，歡迎來信（redraiment@gmail.com）
我的Blogger（子清行）
我的Google Sites（子清行）
我的CSDN博客（夢(mèng)婷軒）
我的百度空間（夢(mèng)婷軒）