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


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,an-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 化為 2s*d 之型式。
Loop K
    隨機產生 [1, n-1] 亂數
    x = ad 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 化為 2s*d 之型式 ,如何將 N-1 化為 2s*d ? 以二進位方式思考很容易想出來。假設
N = 1001001 (N 必為奇數),則 N-1 = 1001000, (N-1)/2 = (N-1) >> 1 = 100100
上式就變成了 100100 = (N-1) = 2s*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;
}

 

這裡示範是只有 100 萬,跑出來結果確實為 78498,而實際上測試另一份數據 - 1億,跑出來結果也確實為 5761445。但再次聲明,這是「機率問題」,將測試次數 TEST_TIMES 調愈高,準確率便愈準。

結束前提醒幾個點:

1. 目前「較新」的質數測試方法為 AKS ,但米勒拉賓測試亦為較常人所知,若對 AKS 有所心得,也請不吝分享討論。

2. 測試質數事實上不該以上述的方式進行使用,因若要找某一區間裡之質數個數,應以篩法為佳。此處只為測試程式之準度而已。

3. 通常質數測試是用來找「大質數」,最小也會超過 100 bits,而不是像這小兒科的小質數。然而若要真正去找「大質數」,勢必會使用 大數演算法。由於本文只撰其概念,故避開了大數演算法之問題。

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

edisonx 發表在 痞客邦 PIXNET 留言(0) 人氣()