[C&++] 深入質數 (1/n) - 試除法單一測試 @ Edison.X. Blog

前二、三個星期幾乎每天都在搞這個，但原始碼實在又多又亂，

說是「深入」其實有點不敢當，畢竟沒有很嚴謹的定義、計算，

只放上一些重點與心得做筆記，也供日後對質數有興趣的網友們一份簡易之參考。

為保持說明之清晰，部份程式碼將不進行任何優化，以最原始方式呈現。

另，若有用到數學符號 log(x)，所代表之意義為 log_e(x)，即為 ln(x)

0. 基本變數宣告

typedef unsigned int uint;
typedef unsigned char byte;

uint i; /* for loop index */
uint n; 
byte IsPrime; /* 當 flag 使用 */

1. 質數之定義

太數學的東西不會說，總之，除了自己與1之外的數，就沒有整數可以整除的話，就是質數。

根據上式定義，可以寫出一份簡明的 code

for(IsPrime = 1, i=2; i<n; ++i) {
   if(n%i==0) {
        IsPrime=0;
        break;
    }
}

但這實在是太慢太慢了。

2. 2n 法則

可以確定的一件事情是，除了2以外，所有的偶數就不再是質數，故在判別時可先濾掉所有的偶數

IsPrime=1;
if(i%2==0) IsPrime=0;
else{
    for(i=3; i<n; i+=2){
       if(n%i==0) {
          IsPrime=0;
          break;
       }
    }
}

這樣約有 50% 左右之數字便不用再判斷

3. 根號法則

有份定理敘述大致是這樣：若一數n ，於根號 n 內仍無其它數字可整除自己，那麼此數便為質數。

有了上面的定理，可再這麼寫下去

IsPrime=1;
if(i%2==0) IsPrime=0;
else{
    for(i=3; i<=(int)ceil(sqrt((double)n)); i+=2){
       if(n%i==0) {
          IsPrime=0;
          break;
       }
    }
}

2n 法則再加上根號法則，想一下，假設要求的數字是 100 萬 (10^6)，實際上要試除的只有 0.5 * sqrt(10^6) = 500 個數字，若再往上，要求1億 (10^8) 以上是不是質數，實際上要試除的數字只有5000 個，1億 <-> 5000，差了二萬倍！

4. 6n 法則

上述說的 2n 法則只有剔除 2 之倍數，要試的是 2n+1，然而 6n 法則是剔除 3 之倍數，實際上只有 6n+1, 6n+5，也就是 6n +- 1 有可能是質數，6個數只要測2個數，這樣下來能省的測試次數從原本 50 % 降到 33%。程式碼大致如下

IsPrime=1;
uint gap=4;
if(n%3==0 || n%2==0) IsPrime=0;
else{
    for(i=5; i<=(int)ceil(sqrt((double)n)); gap=6-gap, i+=gap){
      if(n%i==0) {
          IsPrime=0;
          break;
       }
    }
}

5. 30n 法測

好了，這裡又出現了 30n 法則，這法則我沒看過別人用，後來推斷應是效果也好不到哪去。

2n 法則是剔除 2 之倍數，6n 法則是剔除 2,3 之倍數；30n法則便是剔除 2,3,5這些質數的倍數，

當然也有 210 n 法則，做法都類似。30n 法則該判斷的數為

30n + 1 + 7 +11 +13 +17 +19 +23 +29，

但在實做上沒辦法以 gap= 6-gap 方式去調整，這裡提另一種方式

uint a_s[] = {6,4,2,4,2,4,6,2};//step
uint a_bound = sizeof(a_s)/sizeof(a_s[0])-1;
uint index=0, 
IsPrime=0;
if(n%2==0 || n%3==0 || n%5==0) IsPrime=0;
else{
       for(i=7;  i<=(int)ceil(sqrt((double)n)); i+=a_s[index]){
          if(n%i==0){
             IsPrime=0;
             break;
          }
          index==a_bound ? index=0 : ++index;
       }
}

這種方式明顯比較麻煩，且必須再開一個 array 做 step 動作，

實際上少判斷的數為 8/30 = 26.7%，比 6n 法則少了 6.6 %，

速度也快不了多少，因 array 做存取也需要時間。

6. 210n 法則

而 210n 法則該判斷的數為

    1,  11 ,13,17,19,23,29,31,37,41,43,47,
    53,59,61,67,71,73,79,83,89,97,101,103,107,
    109,113,121,127,131,137,139,143,149,151,157,
    163,167,169,173,179,181,187,191,193,197,199,209

共 48 個；coding 上與上面的 30n 法則具同等效力，於此便不示範 coding。

實際上需判斷數為 48/210 = 22.86%，只比 30n 法則3.84 %，

且實做上也比 30n 法則、6n 法則還慢，目前於二種 CPU 底下測結果，

30n > 6n > 210n。

於此，以 1 億以上之質數為例，以傳統定義方式去求與 30n + 根號法則去比較

一個算 1 億次，一個算 2600 次，試除次數差了三萬八千多倍。

7. 質數表技巧

上述的概念都是以單一試除的概念去判斷一數是否為質數，

若要判斷單一數字是否為質數，可再用一質數表之技巧，步驟大致如下

(1) 程式一開始先開一個較大陣列 (如 10000)，作為質數表。
(2) 先求得一定數量或大小內之質數 (如先求 5000 個質數)。
(3) 日後要判斷某數是否為質數時步驟再如下
      (a) 去質數表裡看該數在不在質數表裡。
      (b) 在做 (a) 時，也做試除動作。
      (c) 步驟 (a), (b) 都無法判斷時，再進行試除
(4) 若求出來為質數，再放到表裡面去。

這技巧只是大致述敘與說明，實做上仍有細節可進行修改與加速。

參考資料

wiki - Mersenne prime
wiki - Fermat number
wiki - Carmichael number

wiki - Sieve of Eratosthenes

wiki - Modular exponentialiation
wiki - Montgomery reduction

wiki - Fermat primality test
wiki - Miller-Rabin primality test
wiki - Lucas Lehmer test for Mersenne numbers
wiki - AKS primality test

wiki - Prime-counting function
wiki - Euler's totient function