[C&++] 深入質數 (3/n) - 質數測試 @ Edison.X. Blog

[警示] 此篇文章所展示之程式碼，只是為「概念上」之程式碼，實際上使用時，會有「溢位」問題。

13. 費馬小定理

質數測試法就我所知主要有幾個：費馬測試(Fermat primality test)、米勒拉賓測試(Miller-Rabin primality test)、盧卡斯測試(Lucas Lehmer test)、AKS test 四種。其中 Lucas 、Fermat 分別判別是否為梅森質數與費馬質數進行之測試，於此不贅述。而 AKS test 需一定數論基礎，這裡跳過 (其實是我也沒數論基礎，看不懂 )，這裡只針對 Fermat primality test / Miller-Rabin primality test 做說明。

費馬小定理，敘述約如下

若 n 為質數，則對所有之互質 a，a^n-1 三 1(mod n)

有幾個邏輯性問題要確認，

(a) n 為質數時， a^(n-1) mod n = 1 必成立。
(b) n 不為質數時，a^(n-1) mod n =1 也可能成立。
(c) 若 a^(n-1) mod n != 1，那 n 就不是質數了。

所以上述是有機率性問題的。若以程式語言表達，大概是長這樣

if ( (uint)pow(a, n-1) % n != 1) printf(" n is not a prime\n");
else printf("n is a prime with prob.\n");

重點便在於，那個 a 要怎麼決定？費馬也給了一段簡單的費馬檢驗，大致如下

for (i=0 ; i!=t; ++i){
a = rand_between(1, n-1);
if ( (uint)pow(a, n) % n != 1) printf(" n is not a prime\n");
}
printf("n is a prime with prob.\n");

注意到，上面的 a 是用隨機去挑，從 1~n-1 去挑的，而最外面的回圈是欲測試之次數測試次數愈多，可信度自然愈高。

然而，有一類型的數會使得費馬測試變得不可靠，此數即為卡麥克數 (Carmichael numbers)。然而卡麥克數非常少，且隨著數字變大，卡麥克數又更少。但能確定的是，卡麥克數都是一些質數之乘積組合而成。但僅管卡麥克數稀少，卻也是使用上，費馬測試不如拉賓測試來得廣泛之原因。

14. 費馬測試

上述對於費馬測試已有初部說明，但實際 coding 起來便又是另一問題，若是要試探一大數是否為質數，該如何測試？以 a^n 為例，若 a=100，n=10000001，這樣用 a^n 必爆，於是必須用一點技巧。

費馬定理裡面提到的是 a^(n-1) mod n = 1，這裡有個 mod 性質可用

(a*b) mod c = ( (a mod c) * ( b mod c) ) mod c

於是 a^n-1 mod n 可簡單的寫下這段 code

unsigned mod_(unsigned a, unsigned n)
{
    unsigned i, r=1;
    for(i=1; i!=n; ++i){
       r = r * a % n;
    }
    return r;
}

然而這段 code 效率並不好，能有更好的方式，與 O(log) 版之 pow 函式相似，此法應稱 Montgomery (蒙格馬利快速取模)

以 p^n % m 為例，程式碼大致如下

unsigned Montgomery(unsigned n, unsigned p, unsigned m)
{    
    unsigned k=1;
    n%=m;
    while(p!=1) {
       if((p&1)!=0) k= (k*n) % m;
       n = (n*n % m);
       p>>=1;
    }
 return (n*k) % m;
}

上面這些都可再替代過，完整程式碼大致如下

#include <stdio.h>
#include <time.h>
#include <stdlib.h>

/* n^p % m */
unsigned Montgomery(unsigned n, unsigned p, unsigned m)
{    
    unsigned k=1;
    n%=m;
    while(p!=1) {
       if((p&1)!=0) k= (k*n) % m;
       n = (n*n % m);
       p>>=1;
    }     
    return (n*k) % m;
}

int fermat_test(unsigned n, unsigned test_times)
{
    
    unsigned i, j, rnd;
// unsigned table[] = { /* 給一定質數表, 這裡列 100 以下 */
// 5, 7, 11, 13, 17, 19, 23, 
// 29, 31, 37, 41, 43, 47, 53, 59, 
// 61, 67, 71, 73, 79, 83, 89, 97 
// };
// const unsigned Size = sizeof(table) / sizeof(table[0]);

    /* 先濾掉 2, 3, 倍數 */
// if(n%2==0 || n%3==0) return 0;

    for(i=0; i!=test_times; ++i){
      rnd = (int)( (double)rand() * (n-4) / RAND_MAX + 2); /* [2,n-2] */
      if(Montgomery(rnd, n-1, n)!=1) return 0; /* 必非質數 */
    }
// /* 可能為質數, 再從質數表裡試除 */
// for(i=0; i!=Size; ++i){
// if(n==table[i]) return 1; /* 直接找此為質數 */
// else if(n%table[i]==0) return 0; /**/
// }
    return 1;
}

int main()
{
    srand((unsigned)time(NULL));
    printf("fermat_test 3215031751:%d\n", fermat_test(3215031751, 1));
    return 0;
}

其中註解部份為保險加上，若為拉賓測試，則註解部份完全都不用理它，

因拉賓測試還有後續動作，有興趣可試試原始版的進行測試。

於此特別註明，也有另一種做法認為，在進行費馬測試時，可直接從質數表裡面隨機挑出質數進行測試，

但原版之費馬測試是聲明 [1, n-1] 隨機挑，故此處從原版 (有不同之處為，從 [2,n-2] 挑 )

15. 拉賓測試

拉賓測試流程還蠻長的，依 wiki 說明大致如下所述

輸入 : N > 3　之奇數，K, 測式之參數，可隨機產生

將 N-1 化為 2^s*d 之型式。
Loop K
    隨機產生 [1, n-1] 亂數
    x = a^d mod n
    if(x==1 || x==n-1) continue;
    for (r=0... s-1)
        x = x*x mod n
        if(x==1) 傳回合數
        if(x==n-1) Next Loop K
    測試完後仍無法判定為質數，傳回合數
Next Loop k
傳回機率質數

重點應在將 N-1 化為 2^s*d 之型式，如何將 N-1 化為 2^s*d ? 以二進位方式思考很容易想出來。假設
N = 1001001 (N 必為奇數)，則 N-1 = 1001000， (N-1)/2 = (N-1) >> 1 = 100100
上式就變成了 100100 = (N-1) = 2^s*d
又根據 shift 原理， 100100 = 1001 << 2，於是 d = 1001， s=2。

有了上述之概念，整份完整程式碼如下

#include <stdio.h>
#include <stdlib.h>
#include <time.h>

#define N 1000
#define TEST_TIMES 10

unsigned random(unsigned low, unsigned up)
{
    unsigned x;
    x = rand();
    // return x % (up-low+1) + low;
    return ((int)( (double)(up-low)*x/RAND_MAX + low ));
}

unsigned Montgomery(unsigned n, unsigned p, unsigned m)
{    
    unsigned k=1;
    n%=m;
    while(p!=1) {
       if((p&1)!=0) k= (k*n) % m;
       n = (n*n % m);
       p>>=1;
    }     
    return (n*k) % m;
}


int MillerTest(unsigned n, unsigned times)
{
    unsigned x, s, d; // N-1 = d*2^s
    unsigned a, r, k;
    unsigned flag=0;

    if(n==2 || n==3) return 1;
    if(!(n&1)) return 0; /* filt multi. of 2 */
    /* n-1 = 2^s*d */
    d=n-1, s=0;
    while(!(d&1)) ++s, d>>=1;
    for(k=0; k!=times; ++k){
       a = random(1, n-1); /* a = rand_between(2, n-2) */
       x = Montgomery(a, d, n); /* x=a^d mod n */
       if(x==1 || x==n-1) continue;
       for(r=0; r!=s; ++r){
          if(x==n-1) break;
          // x = Montgomery(x, 2, n);
          x = (x*x) % n;
       }
       if(r<s) continue;
       return 0;
    }
    return 1;
}

int main()
{
    unsigned i, cnt=0;

    srand((unsigned)time(NULL));
    for(i=2; i!=N; ++i) cnt+=MillerTest(i, TEST_TIMES);
    printf("cnt=%d\n", cnt);
    return 0;
}

好了，拿去測試之後，N=10, 100, 1000, 10000 都很準。疑，那 10 萬呢？

問題就在這裡，10 萬時，出來的數字反而開始不準。10 萬以下之亂數應有 9592 個，但實際上測得之亂數個數約只有 6579 個，為何會這樣？

原因應有二個 - 亂數產生器、overflow。我手邊之亂數產生器最大只到 32767，若要再上去，就必須用 ( rand() << 15) | rand() 方式，若怕不均勻的話，可再去找其它亂數產生器。而另一問題應出在取模函式裡面，因裡面用到了 a*m %n 這東西，若在取模之前便 overflow，那取模式沒意義。這部份倒是可用 (a*m) % n = ( (a%n) * (m%n) ) % n; 加以修正，但這麼做無庸便多用了二次取模，效能變差。

另還有個簡單方式可解決，但要附出一些相對成本，即為將資料型態改成 unsigned long long ，此時在 random 便不適用 return ((uint64)( (double)(up-low)*x/RAND_MAX + low )); 此方式，因 double 本身有精度問題，無法完全表示 unsigned long long 範圍，故仍以基本之取模運算取得。

另一需注意，我手邊電腦 RAND_MAX 為 32767，15 bits，要產生夠大的亂數，可考慮這麼做

uint64 x = ( (uint64)rand() & 0xFFF0000000000000ULL >> 0 ) |
                 ( (uint64)rand() & 0xFFF1000000000000ULL >> 12 ) |
                 ( (uint64)rand() & 0xFFF1000000000000ULL >> 25 ) |
                 ( (uint64)rand() & 0xFFF1000000000000ULL >> 38 ) |
                 ( (uint64)rand() & 0xFFF1000000000000ULL >> 41 ) ;

上述程式碼看起來麻煩些，但就亂數品質而言，是較下面方式佳一點，避免使用以下方式 (感謝 Jacob 給予指正與意見)

uint64 x = ( rand() << 45 ) | ( rand() << 30 ) | ( rand() << 15 ) | rand(); /* worst */

這部份有機會再開「亂數」文章予以探討。

如上所言，這麼做的話這隻亂數產生器品質已經變壞很多了，最好的方式還是再生一隻亂數產生器出來，如 C++ 裡可調用 std::tr1 使用。修改過後程式碼如下所示

#include <stdio.h>
#include <stdlib.h>
#include <time.h>

#define N 100000
#define TEST_TIMES 10

typedef unsigned long long uint64;

uint64 random(uint64 low, uint64 up)
{
    uint64 x = ( (uint64)rand() & 0xFFF0000000000000ULL >> 0 ) |
       ( (uint64)rand() & 0xFFF1000000000000ULL >> 12 ) |
       ( (uint64)rand() & 0xFFF1000000000000ULL >> 25 ) |
       ( (uint64)rand() & 0xFFF1000000000000ULL >> 38 ) |
       ( (uint64)rand() & 0xFFF1000000000000ULL >> 41 ) ;

    return x % (up-low+1) + low;
}

uint64 Montgomery(uint64 n, uint64 p, uint64 m)
{    
    uint64 k=1;
    n%=m;
    while(p!=1) {
       if((p&1)!=0) k= (k*n) % m;
       n = (n*n % m);
       p>>=1;
    }     
    return (n*k) % m;
}


int MillerTest(uint64 n, uint64 times)
{
    uint64 x, s, d; // N-1 = d*2^s
    uint64 a, r, k;
    uint64 flag=0;

    if(n==2 || n==3) return 1;
    if(!(n&1)) return 0; /* filt multi. of 2 */
    /* n-1 = 2^s*d */
    d=n-1, s=0;
    while(!(d&1)) ++s, d>>=1;
    for(k=0; k!=times; ++k){
       a = random(1, n-1); /* a = rand_between(2, n-2) */
       x = Montgomery(a, d, n); /* x=a^d mod n */
       if(x==1 || x==n-1) continue;
       for(r=0; r!=s; ++r){
          if(x==n-1) break;
          // x = Montgomery(x, 2, n);
          x = (x*x) % n;
       }
       if(r<s) continue;
       return 0;
    }
    return 1;
}

int main()
{
    uint64 i;
    int cnt=0;

    srand((uint64)time(NULL));
    for(i=2; i!=N; ++i) cnt+=MillerTest(i, TEST_TIMES);
    printf("cnt=%d\n", cnt);
    return 0;
}