深入淺出素數演算法

// 缺少的代碼請參考前邊</p><p>#include <stdlib.h></p><p>static int cmp(const int *p, const int *q)<br />{<br /> return (*p) - (*q);<br />}</p><p>bool isPrime(int n)<br />{<br /> if(n < 2) return false;<br /> if(n == 2) return true;<br /> if(n%2 == 0) return false;</p><p> if(n >= 67 && n <= primes[NELEMS(primes)-1])<br /> {<br /> return NULL !=<br /> bsearch(&n, primes, NELEMS(primes), sizeof(n), cmp);<br /> }<br /> else<br /> {<br /> for(int i = 1; primes[i]*primes[i] <= n; ++i)<br /> if(n%primes[i] == 0) return false;<br /> return true;<br /> }<br />}<br />

注意: 如果沒有特殊說明, 以下討論的都是針對n為素數時的時間複雜度

1. 根據概念判斷:

     如果一個正整數只有兩個因子, 1和p，則稱p為素數.

bool isPrime(int n)<br />{<br /> if(n < 2) return false;</p><p> for(int i = 2; i < n; ++i)<br /> if(n%i == 0) return false;</p><p> return true;<br />}</p><p>

時間複雜度O(n).

2. 改進, 去掉偶數的判斷

bool isPrime(int n)<br />{<br /> if(n < 2) return false;<br /> if(n == 2) return true;</p><p> for(int i = 3; i < n; i += 2)<br /> if(n%i == 0) return false;</p><p> return true;<br />}</p><p>

時間複雜度O(n/2), 速度提高一倍.

3. 進一步減少判斷的範圍

     定理: 如果n不是素數, 則n有滿足1<d<=sqrt(n)的一個因子d.
     證明: 如果n不是素數, 則由定義n有一個因子d滿足1<d<n.
             如果d大於sqrt(n), 則n/d是滿足1<n/d<=sqrt(n)的一個因子.

bool isPrime(int n)<br />{<br /> if(n < 2) return false;<br /> if(n == 2) return true;</p><p> for(int i = 3; i*i <= n; i += 2)<br /> if(n%i == 0) return false;</p><p> return true;<br />}</p><p>

這樣也可以

bool IsPrim(int a)<br />{<br /> int divisor=3;<br /> int limit=a;<br /> if(a%2==0)<br /> return false;<br /> while(limit>divisor)<br /> {<br /> if(a%divisor==0)<br /> return false;<br /> limit=a/divisor;<br /> divisor+=2;<br /> }<br /> return true;<br />}                               
      
  

時間複雜度O(sqrt(n)/2), 速度提高O((n-sqrt(n))/2).

4. 剔除因子中的重複判斷.
     例如: 11%3 != 0 可以確定 11%(3*i) != 0.

     定理: 如果n不是素數, 則n有滿足1<d<=sqrt(n)的一個"素數"因子d.
     證明: I1. 如果n不是素數, 則n有滿足1<d<=sqrt(n)的一個因子d.
             I2. 如果d是素數, 則定理得證, 演算法終止.
             I3. 令n=d, 並轉到步驟I1.

     由於不可能無限分解n的因子, 因此上述證明的演算法最終會停止.

// primes[i]是遞增的素數序列: 2, 3, 5, 7, ...<br />// 更準確地說primes[i]序列包含1->sqrt(n)範圍內的所有素數</p><p>bool isPrime(int primes[], int n)<br />{<br /> if(n < 2) return false;</p><p> for(int i = 0; primes[i]*primes[i] <= n; ++i)<br /> if(n%primes[i] == 0) return false;</p><p> return true;<br />}</p><p>

假設n範圍內的素數個數為PI(n), 則時間複雜度O(PI(sqrt(n))).

函數PI(x)滿足素數定理: ln(x)-3/2 < x/PI(x) < ln(x)-1/2, 當x >= 67時.

因此O(PI(sqrt(n)))可以表示為O(sqrt(x)/(ln(sqrt(x))-3/2)),

O(sqrt(x)/(ln(sqrt(x))-3/2))也是這個演算法的空間複雜度.

5. 構造素數序列primes[i]: 2, 3, 5, 7, ...

由4的演算法我們知道, 在素數序列已經被構造的情況下, 判斷n是否為素數效率很高;

但是, 在構造素數序列本身的時候, 是否也可是達到最好的效率呢?

事實上這是可以的! -- 我們在構造的時候完全可以利用已經被構造的素數序列!

假設我們已經我素數序列: p1, p2, .. pn

現在要判斷pn+1是否是素數, 則需要(1, sqrt(pn+1)]範圍內的所有素數序列,

而這個素數序列顯然已經作為p1, p2, .. pn的一個子集被包含了!

// 構造素數序列primes[]</p><p>void makePrimes(int primes[], int num)<br />{<br /> int i, j, cnt;</p><p> primes[0] = 2;<br /> primes[1] = 3;</p><p> for(i = 5, cnt = 2; cnt < num; i += 2)<br /> {<br /> int flag = true;<br /> for(j = 1; primes[j]*primes[j] <= i; ++j)<br /> {<br /> if(i%primes[j] == 0)<br /> {<br /> flag = false; break;<br /> }<br /> }<br /> if(flag) primes[cnt++] = i;<br /> }<br />}</p><p>

makePrimes的時間複雜度比較複雜, 而且它只有在初始化的時候才被調用一次.在一定的應用範圍內, 我們可以把近似認為makePrimes需要常數時間.在後面的討論中, 我們將探討一種對電腦而言更好的makePrimes方法.

附：素數的刪法

     [定理]若比素數P小的所有素數的倍數均已從IsPrim中刪去，且P*P > N,則 剩下的數就全為素數。

 void CreatPrim()<br /> {<br /> for(int i=2;i<=N;i++)<br /> IsPrim[i]=i%2;<br /> IsPrim[2]=1;<br /> for(int i=3;i*i<N;i+=2)<br /> {<br /> if(IsPrim[i])<br /> for(int j=i;i*j<=N;j+=2)<br /> IsPrim[i*j]=0;<br /> }<br /> } </p><p>
6. 更好地利用電腦資源...

當前的主流PC中, 一個整數的大小為2^32. 如果需要判斷2^32大小的數是否為素數,則可能需要測試[2, 2^16]範圍內的所有素數(2^16 == sqrt(2^32)).由4中提到的素數定理我們可以大概確定[2, 2^16]範圍內的素數個數.由於2^16/(ln(2^16)-1/2) = 6138, 2^16/(ln(2^16)-3/2) = 6834,我們可以大概估計出[2, 2^16]範圍內的素數個數6138 < PI(2^16) < 6834.

在對[2, 2^16]範圍內的素數進行統計, 發現只有6542個素數:

p_6542: 65521, 65521^2 = 4293001441 < 2^32, (2^32 = 4294967296)
p_6543: 65537, 65537^2 = 4295098369 > 2^32, (2^32 = 4294967296)

在實際運算時unsigned long x = 4295098369;將發生溢出, 為131073.在程式中, 我是採用double類型計算得到的結果.

分析到這裡我們可以看到, 我們只需要緩衝6543個素數, 我們就可以採用4中的演算法高效率地判斷[2, 2^32]如此龐大範圍內的素數!(原本的2^32大小的問題規模現在已經被減小到6543規模了!)

雖然用現在的電腦處理[2, 2^16]範圍內的6542個素數已經沒有一點問題,雖然makePrimes只要被運行一次就可以, 但是我們還是考慮一下是否被改進的可能?!

我想學過java的人肯定想把makePrimes作為一個靜態初始化實現, 在C++中也可以類比java中靜態初始化的類似實現:

#define NELEMS(x) ((sizeof(x)) / (sizeof((x)[0])))

static int primes[6542+1];
static struct _Init { _Init(){makePrimes(primes, NELEMS(primes);} } _init;

如此, 就可以在程式啟動的時候自動掉用makePrimes初始化素數序列.但, 我現在的想法是: 為什麼我們不能在編譯的時候調用makePrimes函數呢? 完全可以!!! 代碼如下:

// 這段代碼可以由程式直接產生</p><p>const static int primes[] =<br />{<br />2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,101,103,<br />107,109,113,127,131,137,139,149,151,157,163,167,173,179,181,191,193,197,199,211,<br />223,227,229,233,239,241,251,257,263,269,271,277,281,283,293,307,311,313,317,331,<br />337,347,349,353,359,367,373,379,383,389,397,401,409,419,421,431,433,439,443,449,<br />457,461,463,467,479,487,491,499,503,509,521,523,541,547,557,563,569,571,577,587,<br />593,599,601,607,613,617,619,631,641,643,647,653,659,661,673,677,683,691,701,709,<br />719,727,733,739,743,751,757,761,769,773,787,797,809,811,821,823,827,829,839,853,<br />857,859,863,877,881,883,887,907,911,919,929,937,941,947,953,967,971,977,983,991,<br />...<br />65521, 65537<br />};</p><p>

有點不可思議吧:), 原本makePrimes需要花費的時間複雜度現在真的變成O(1)了!(我覺得叫O(0)可能更合適!)

7. 二分法尋找

現在我們緩衝了前大約sqrt(2^32)/(ln(sqrt(2^32)-3/2))個素數列表, 在判斷2^32層級的

素數時最多也只需要PI(sqrt(2^32))次判斷(準確值是6543次), 但是否還有其他的方式判斷呢?

當素數比較小的時候(不大於2^16), 是否可以直接從緩衝的素數列表中直接查詢得到呢?

答案是肯定的! 由於primes是一個有序的數列, 因此我們當素數小於2^16時, 我們可以直接

採用二分法從primes中查詢得到(如果查詢失敗則不是素數).

 

時間複雜度:

    if(n <= primes[NELEMS(primes)-1] && n >= 67): O(log2(NELEMS(primes))) < 13;
    if(n >    primes[NELEMS(primes)-1]): O(PI(sqrt(n))) <= NELEMS(primes).

8. 素數定理+2分法尋找

在9中, 我們對小等於primes[NELEMS(primes)-1]的數採用2分法尋找進行判斷.我們之前針對2^32緩衝的6453個素數需要判斷的次數為 13次(log2(1024*8) == 13).對於小的素數而言(其實就是2^16範圍只內的數), 13次的比較已經完全可以接受了.不過根據素數定理: ln(x)-3/2 < x/PI(x) < ln(x)-1/2, 當x >= 67時, 我們依然可以進不步縮小小於2^32情況的尋找範圍(現在是0到NELEMS(primes)-1範圍尋找).我們需要解決問題是n <= primes[NELEMS(primes)-1):

如果n為素數, 那麼它在素數序列可能出現的範圍在哪?

      ---- (n/(ln(n)-1/2), n/(ln(n)-3/2)), 即素數定理!

上面的代碼修改如下:

bool isPrime(int n)<br />{<br /> if(n < 2) return false;<br /> if(n == 2) return true;<br /> if(n%2 == 0) return false;</p><p> int hi = (int)ceil(n/(ln(n)-3/2));</p><p> if(n >= 67 && hi < NELEMS(primes))<br /> {<br /> int lo = (int)floor(n/(ln(n)-1/2));</p><p> return NULL !=<br /> bsearch(&n, primes+lo, hi-lo, sizeof(n), cmp);<br /> }<br /> else<br /> {<br /> for(int i = 1; primes[i]*primes[i] <= n; ++i)<br /> if(n%primes[i] == 0) return false;<br /> return true;<br /> }<br />}</p><p>

時間複雜度:

    if(n <= primes[NELEMS(primes)-1] && n >= 67): O(log2(hi-lo))) < ???;
    if(n >    primes[NELEMS(primes)-1]): O(PI(sqrt(n))) <= NELEMS(primes).

9. 打包成素數庫(給出全部的代碼)

到目前為止, 我已經給出了我所知道所有改進的方法(如果有人有更好的演算法感謝告訴我).這裡需要強調的一點是, 這裡討論的素數求法是針對0-2^32範圍的數而言, 至於像尋找成百上千位大小的數不在此討論範圍, 那應該算是純數學的內容了.

代碼儲存在2個檔案: prime.h, prime.cpp.

// file: prime.h</p><p>#ifndef PRIME_H_2006_10_27_<br />#define PRIME_H_2006_10_27_</p><p>extern int Prime_max(void); // 素數序列的大小<br />extern int Prime_get (int i); // 返回第i個素數, 0 <= i < Prime_max</p><p>extern bool Prime_test(int n); // 測試是否是素數, 1 <= n < INT_MAX</p><p>#endif</p><p>///////////////////////////////////////////////////////</p><p>// file: prime.cpp</p><p>#include <assert.h><br />#include <limits.h><br />#include <math.h><br />#include <stdlib.h></p><p>#include "prime.h"</p><p>// 計算數組的元素個數</p><p>#define NELEMS(x) ((sizeof(x)) / (sizeof((x)[0])))</p><p>// 素數序列, 至少儲存前6543個素數!<br />// 我個人習慣緩衝前（1024*8）個:-)</p><p>static const int primes[] =<br />{<br />2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,101,103,<br />107,109,113,127,131,137,139,149,151,157,163,167,173,179,181,191,193,197,199,211,<br />223,227,229,233,239,241,251,257,263,269,271,277,281,283,293,307,311,313,317,331,<br />337,347,349,353,359,367,373,379,383,389,397,401,409,419,421,431,433,439,443,449,<br />457,461,463,467,479,487,491,499,503,509,521,523,541,547,557,563,569,571,577,587,<br />593,599,601,607,613,617,619,631,641,643,647,653,659,661,673,677,683,691,701,709,<br />719,727,733,739,743,751,757,761,769,773,787,797,809,811,821,823,827,829,839,853,<br />857,859,863,877,881,883,887,907,911,919,929,937,941,947,953,967,971,977,983,991,<br />...<br />65521, 65537<br />};</p><p>// bsearch的比較函數</p><p>static int cmp(const void *p, const void *q)<br />{<br /> return (*(int*)p) - (*(int*)q);<br />}</p><p>// 緩衝的素數個數</p><p>int Prime_max()<br />{<br /> return NELEMS(primes);<br />}</p><p>// 返回第i個素數</p><p>int Prime_get(int i)<br />{<br /> assert(i >= 0 && i < NELEMS(primes));<br /> return primes[i];<br />}</p><p>// 測試n是否是素數</p><p>bool Prime_test(int n)<br />{<br /> assert(n > 0);</p><p> // 偶數情況單獨判斷</p><p> if(n < 2) return false;<br /> if(n == 2) return true;<br /> if(!(n&1)) return false;</p><p> // 如果n為素數, 則在序列hi位置之前</p><p> int lo, hi = (int)ceil(n/(log(n)-3/2.0));</p><p> if(hi < NELEMS(primes))<br /> {<br /> // 確定2分法尋找的範圍<br /> // 只有n >= 67是才滿足素數定理</p><p> if(n >= 67) lo = (int)floor(n/(log(n)-1/2.0));<br /> else { lo = 0; hi = 19; }</p><p> // 尋找成功則為素數</p><p> return NULL !=<br /> bsearch(&n, primes+lo, hi-lo, sizeof(n), cmp);<br /> }<br /> else<br /> {<br /> // 不在儲存的素數序列範圍之內的情況</p><p> for(int i = 1; primes[i]*primes[i] <= n; ++i)<br /> if(n%primes[i] == 0) return false;</p><p> return true;<br /> }<br />}</p><p>
10. 回顧, 以及推廣

到這裡, 關於素數的討論基本告一段落. 回顧我們之前的求解過程, 我們會發現如果缺少數學的基本知識會很難設計好的演算法; 但是如果一味地只考慮數學原理,而忽律了電腦的本質特徵, 也會有同樣的問題.一個很常見的例子就是求Fibonacci數列. 當然方法很多, 但是在目前的電腦中都沒有實現的必要!

因為Fibonacci數列本身是指數增長的, 32位的有符號整數所能表示的位置只有前46個:

static const int Fibonacci[] =<br />{<br /> 0,1,1,2,3,5,8,13,21,34,55,89,144,233,377,610,987,1597,<br /> 2584,4181,6765,10946,17711,28657,46368,75025,121393,196418,<br /> 317811,514229,832040,1346269,2178309,3524578,5702887,9227465,<br /> 14930352,24157817,39088169,63245986,102334155,165580141,267914296,<br /> 433494437,701408733,1134903170,1836311903,-1323752223</p><p> // 注意：數列到F47已經發生溢出！！！<br />};</p><p>

因此, 我只需要把前46個Fibonacci數儲存到數組中就可以搞定了!比如: F(int i){return Fibonacci(i);}非常簡單, 效率也非常好.同樣的例子如求階乘n!, 雖然也有很多數學上的描述, 但是電腦一般都表示不了,我們直接把能用的計算好放到記憶體中就可以了.

總之, 許多東西本身是好的, 但是不要被它束縛了!