斐波那契序列 集錦 (轉)
[定理1] 標準Fibonacci序列（即第0項為0，第1項為1的序列）當N大於1時，一定有f(N)和f(N-1)互質

其實，結合“互質”的定義，和一個很經典的演算法就可以輕鬆證明 
對，就是輾轉相除法 
互質的定義就是最大公約數為1

數學歸納法是很有用的證明方法，我們接下來這個定理用數學歸納法就很好證明： 
[定理2]若i為奇數， f(i)*f(i)=f(i-1)*f(i+1)+1，否則f(i)*f(i)=f(i-1)*f(i+1)-1 
對，這個定理用數學歸納法可以輕鬆證明，大家有興趣可以自己嘗試

[定理3] f(n)=f(i)*f(n-i-1)+f(i+1)*f(n-i)

f(n)=f(1)*f(n-2)+ f(2)*f(n-1) 
    =f(2)*f(n-3)+ f(3)*f(n-2) 
    =f(3)*f(n-4)+ f(4)*f(n-3) 
看來沒有錯，證明方法就是這樣

這個公式也可以用來計算較大的fibonacci數除以某個數的餘數

設i=n/2 不過，為了保證計算能延續下去 需要每次保留三個值 
這樣，下一次計算就可以利用這三個值來求出兩個值，再相加就可以得到第三個值 
譬如，計算出f(5),f(6),f(7)，可以計算出f(11)、f(12)，然後推出f(13) 
就是剛才洛奇的悲鳴(364738334)所提到的矩陣方法 
我們知道我們若要簡單計算f(n)，有一種方法就是先儲存 
a=f(0),b=f(1)，然後每次設： 
a'=b b'=a+b

並用新的a'和b'來繼續這一運算

如果大家熟悉利用“矩陣”這一工具的話，就知道，如果把a、b寫成一個向量[a,b]，完成上述操作相當於乘以矩陣 
0 1 
1 1 
也就是說，如果我們要求第100個fibonacci數，只需要將矩陣 
[0,1]乘上 
0 1 
1 1 
的一百次方，再取出第一項

因為我們知道，矩陣運算滿足結合律，一次次右乘那個矩陣完全可以用乘上那個矩陣的N次方代替，更進一步，那個矩陣的N次方就是這樣的形式： 
f(n-1) f(n) 
f(n) f(n+1)

而求矩陣的N次方，由於矩陣乘法滿足結合律，所以我們可以用log(N)的演算法求出——這個演算法大家都會麼？ 
一個是二分，一個是基於二進位制的求冪

二分的原理：要求矩陣的N次方A(N)，設i=N/2若N%2==1， 則 A(N)=A(i)*A(i)*A(1)若N%2==0， 則 A(N)=A(i)*A(i)

基於二進位制的原理：將N拆為二進位制數，譬如13=1101那麼 A^13= A^8 * A^4 * A^1 （這裡^表示冪運算）

也就是說，由A^1開始，自乘得到A^2，然後自乘得到A^4，如果N對應位為1，則將這個結果乘到目標上去

這樣的話，將所有乘法改為模乘，就可以得到一個較大Fibonacci數除以M的餘數

若不用遞迴，其實類似

這題是要求求F(n)的最後四位數，所有乘法過程增加一個模10000的步驟即可，大家可以收藏稍候AC

關於矩陣我們告一段落，等下會回來繼續探討利用矩陣來解決複雜些的Fibonacci問題

當然，用矩陣也可以解決這道題——只要將乘法改為乘並保留前四位

我們採用double 保留整數部分四位 這題最好還是double吧

不過顯然有更好的解法——如果我們知道Fibonacci序列的通項公式

F(n) = (((1+Sqrt(5))/2)^n - ((1-Sqrt(5))/2)^n)*1/Sqrt(5)

不過組合數學裡也有這一公式的推導方法 叫做“線性齊次遞推式”

這個解法的核心是，通解是某個數的冪 將f(n)=x^n代入遞推方程，可以解出三個通解 0和 (1+sqrt(5))/2

通常把“0”稱作平凡解，那麼特解就是通解的某個線性組合

再代入f(0)=0 f(1)=1，就可以得出我們剛才的公式

不過通常情況下，我們只需要記住那個公式就可以了

提醒大家，記憶公式的時候千萬別忘記了係數1/sqrt(5)

因為(1-sqrt(5))/2的絕對值小於1

所以當i較大的時候，往往可以忽略掉這一項 
f(i)≈((1+Sqrt(5))/2)^n/sqrt(5);

所以，剛才列舉出的HDOJ的1568，可以很簡單的30以內直接求解，30以上採用這個公式，還是用log(N)求冪的演算法求解 
恩，就是公式的前半部分

[定理5] 標準Fibonacci序列對任意大於2的正整數的餘數序列，必然是以“0 1”為迴圈節開頭的序列

顯然0、1是序列開頭，也就是說序列開頭就是迴圈節開頭

迴圈長度的計算貌似是個比較難的問題，我一時還沒有想到有效解法，不過，要說明的是，計算複雜度時，這個迴圈節長度應該按複雜度O(N^2)計算

恩，證明方法是利用同餘定理、反證法，還有我們之前證明過的相鄰項一定互質的定理，留給大家家庭作業

現在告訴大家一種正確解法，然後大家就可以去搞定這道題向別人炫耀了

首先，我們將問題整理一下，就是對等差數列 ai=k*i+b，求所有的f(ai)之和除以M的餘數

當0<=i<N

大家有沒有想到，因為ai是等差數列，倘若f(ai)也是個等什麼序列，那說不定就有公式求了

f(ai)顯然不是等差數列，直接看上去也不是等比數列

但是如果把f(ai)換成我們剛才所說的Fibonacci矩陣呢？

是的，可是我們對矩陣是直接求冪即可，由於矩陣加法的性質，我們要求A^ai的右上角元素之和，只要求A^ai之和的右上角元素

就矩陣這個東西來說，完全可以看作一個等比數列， 
首項是：A^b 
公比是：A^k 
項數是：N

呵呵，我們可以把問題進一步簡化

因為矩陣的加法對乘法也符合分配律，我們提出一個A^b來，形成這樣的式子： 
A^b*( I + A^k + (A^k)^2 + .... + (A^k)^(N-1) )

A^b 和 A^k 顯然都可以用我們之前說過的方法計算出來，這剩下一部分累加怎麼解決呢

簡單起見，設A^k=B 
要求 G(N)=I + ... + B^(N-1)，設i=N/2 
若N為偶數，G(N)=G(i)+G(i)*B^i 
若N為奇數，G(N)=I+ G(i)*B + G(i) * (B^(i+1))

呵呵，這個方法就是比賽當時ACRush用的 
而農夫用的則是更精妙的方法，大家可想知道

我們來設定這樣一個矩陣 
B I 
O I 
其中O是零矩陣，I是單位矩陣

將它乘方，得到 
B^2 I+B 
O   I 
乘三方，得到 
B^3 I+B+B^2 
O   I 
乘四方，得到 
B^4 I+B+B^2+B^3 
O   I

既然已經轉換成矩陣的冪了，繼續用我們的二分或者二進位制法，直接求出冪就可以了

Fibinary數是指沒有相鄰的兩個1的二進位制數。給N，求出第N大的Fibinary數

相對於二進位制中每一位的值是2的冪，十進位制中每一位的值是十的冪， 
Fibonacci進位制是每一位的值是對應Fibonacci數的一種計數系統。 
     8 5 3 2 1 
1     1 
2     1 0 
3     1 0 0 
4     1 0 1 
5     1 0 0 0 
6     1 0 0 1 
7     1 0 1 0 
8     1 0 0 0 0 
9     1 0 0 0 1 
10   1 0 0 1 0 
11   1 0 1 0 0 
12   1 0 1 0 1 
以上是前12個數字對應的十進位制到Fibonacci進位制的表格

Fibonacci的運算方法很奇怪。首先，它每一位上非0即1，而且不同於二進位制的逢二進一或者十進位制的逢十進一，它的進位方法是逢連續兩個1，則進1

譬如 
1010110==> 1011000 ==> 1100000==>10000000

在最低位有個特殊情況，最低位既可以逢2進1，也可以和次低位一起逢相鄰進1 
這種奇怪的進位方法，換句話描述就是，不存在兩個連續的1 
因為Fibonacci數其實也增長很快，int範圍內好像只有46個，本題只需要用最簡單的辦法轉換成Fibonacii進位制即可 
其中一題是 
http://www.mydrs.org/program/down/ ahoi2004da y1.pdf 
中的第二題，叫做數字迷陣 
還有一題是PKU上的很出名 的取石子問題 
http://acm.pku.edu.cn/JudgeOnline/problem?id=1067

另外這題 可以利用Fibonacci判斷資料範圍進行優化設計。不過貌似可以水過去，僅僅給大家提供個思路罷

posted @ 2009-09-04 00:09 Knuth_檔案 閱讀(50) | 評論 (0) | 編輯

費馬小定理 素數判定 蒙哥馬利演算法 
約定：
x%y為x取模y，即x除以y所得的餘數，當x<y時，x%y=x，所有取模的運算物件都為整數。
x^y表示x的y次方。
乘方運算的優先順序高於乘除和取模，加減的優先順序最低。
見到x^y/z這樣，就先算乘方，再算除法。
A/B，稱為A除以B，也稱為B除A。
若A%B=0，即稱為A可以被B整除，也稱B可以整除A。
A*B表示A乘以B或稱A乘B，B乘A，B乘以A……都TMD的一樣，靠！

複習一下小學數學
公因數：兩個不同的自然數A和B，若有自然數C可以整除A也可以整除B，那麼C就是A和B的公因數。
公倍數：兩個不同的自然數A和B，若有自然數C可以被A整除也可以被B整除，那麼C就是A和B的公倍數。
互質數：兩個不同的自然數，它們只有一個公因數1，則稱它們互質。

費馬是法國數學家，又譯“費爾馬”，此人巨牛，他的簡介請看下面。不看不知道，一看嚇一跳。
　　

費馬小定理：
有N為任意正整數，P為素數，且N不能被P整除（顯然N和P互質），則有：
N^P%P=N(即：N的P次方除以P的餘數是N)

但是我查了很多資料見到的公式都是這個樣子：
(N^(P-1))%P=1

後來分析了一下，兩個式子其實是一樣的，可以互相變形得到，原式可化為：
(N^P-N)%P=0(即：N的P次方減N可以被P整除，因為由費馬小定理知道N的P次方除以P的餘數是N)

把N提出來一個，N^P就成了你N*(N^(P-1))，那麼(N^P-N)%P=0可化為：(N*(N^(P-1)-1))%P=0
請注意上式，含義是：N*(N^(P-1)-1)可以被P整除

又因為N*(N^(P-1)-1)必能整除N（這不費話麼!）
所以，N*(N^(P-1)-1)是N和P的公倍數，小學知識了^_^

又因為前提是N與P互質，而互質數的最小公倍數為它們的乘積，所以一定存在正整數M使得等式成立：
N*(N^(P-1)-1)=M*N*P
兩邊約去N，化簡之：
N^(P-1)-1=M*P
因為M是整數，顯然：
(N^(P-1)-1)%P=0
即：
N^(P-1)%P=1
============================================
積模分解公式

先有一個引理，如果有：X%Z=0，即X能被Z整除，則有：
(X+Y)%Z=Y%Z
這個不用證了吧...

設有X、Y和Z三個正整數，則必有：(X*Y)%Z=((X%Z)*(Y%Z))%Z

想了很長時間才證出來，要分情況討論才行：

1.當X和Y都比Z大時，必有整數A和B使下面的等式成立：
X=Z*I+A（1）
Y=Z*J+B（2）
不用多說了吧，這是除模運算的性質！
將（1）和（2）代入(X*Y)modZ得：((Z*I+A)(Z*J+B))%Z
乘開，再把前三項的Z提一個出來，變形為：(Z*(Z*I*J+I*A+I*B)+A*B)%Z（3）
因為Z*(Z*I*J+I*A+I*B)是Z的整數倍……暈，又來了。
概據引理，（3）式可化簡為：(A*B)%Z
又因為：A=X%Z，B=Y%Z，代入上面的式子，就成了原式了。

2.當X比Z大而Y比Z小時，一樣的轉化：
X=Z*I+A
代入(X*Y)%Z得：
(Z*I*Y+A*Y)%Z
根據引理，轉化得：(A*Y)%Z
因為A=X%Z，又因為Y=Y%Z，代入上式，即得到原式。
同理，當X比Z小而Y比Z大時，原式也成立。

3.當X比Z小，且Y也比Z小時，X=X%Z，Y=Y%Z，所以原式成立。
=====================================================
快速計算乘方的演算法

如計算2^13，則傳統做法需要進行12次乘法。

/*計算n^p*/
unsigned power(unsigned n,unsigned p)
{
    for(int i=0;i<p;i++) n*=n;
    return n;
}

該死的乘法，是時候優化一下了！把2*2的結果儲存起來看看，是不是成了：4*4*4*4*4*4*2 
再把4*4的結果儲存起來：16*16*16*2 
一共5次運算，分別是2*2、4*4和16*16*16*2

這樣分析，我們演算法因該是隻需要計算一半都不到的乘法了。
為了講清這個演算法，再舉一個例子2^7：2*2*2*2*2*2*2 
兩兩分開：(2*2)*(2*2)*(2*2)*2 
如果用2*2來計算，那麼指數就可以除以2了，不過剩了一個，稍後再單獨乘上它。
再次兩兩分開，指數除以2： ((2*2)*(2*2))*(2*2)*2 
實際上最後一個括號裡的2 * 2是這回又剩下的，那麼，稍後再單獨乘上它 
現在指數已經為1了，可以計算最終結果了：16*4*2=128

優化後的演算法如下：
unsigned Power(unsigned n,unsigned p) 
{
   unsigned main=n; //用main儲存結果
   unsigned odd=1; //odd用來計算“剩下的”乘積
   while (p>1) 
   {//一直計算，直到指數小於或等於1
        if((p%2)!=0)
        {// 如果指數p是奇數，則說明計算後會剩一個多餘的數，那麼在這裡把它乘到結果中
            odd*=main; //把“剩下的”乘起來
        }
        main*=main; //主體乘方
        p/=2; //指數除以2
   }
   return main*odd; //最後把主體和“剩下的”乘起來作為結果
}

夠完美了嗎？不，還不夠！看出來了嗎？main是沒有必要的，並且我們可以有更快的程式碼來判斷奇數。要知道除法或取模運算的效率很低，所以我們可以利用偶數的一個性質來優化程式碼，那就是偶數的二進位制表示法中的最低位一定為0!

完美版：
unsigned Power(unsigned n, unsigned p) 
{ // 計算n的p次方
    unsigned odd = 1; //oddk用來計算“剩下的”乘積
    while (p > 1)
    { // 一直計算到指數小於或等於1
       if (( p & 1 )!=0)
      { // 判斷p是否奇數，偶數的最低位必為0
             odd *= n; // 若odd為奇數，則把“剩下的”乘起來
      }
      n *= n; // 主體乘方
      p /= 2; // 指數除以2
     }
    return n * odd; // 最後把主體和“剩下的”乘起來作為結果
}
========================================================
蒙格馬利”快速冪模演算法

後面我們會用到這樣一種運算：(X^Y)%Z

問題是當X和Y很大時，只有32位的整型變數如何能夠有效的計算出結果？
考慮上面那份最終的優化程式碼和再上面提到過的積模分解公式，我想你也許會猛拍一下腦門，吸口氣說：“哦，我懂了！”。

下面的講解是給尚沒有做出這樣動作的同學們準備的。X^Y可以看作Y個X相乘，即然有積模分解公式，那麼我們就可以把Y個X相乘再取模的過程分解開來，比 如：(17^25)%29則可分解為：( ( 17 * 17 ) % 29 * ( 17 * 17 ) % 29 * ……
如果用上面的程式碼將這個過程優化，那麼我們就得到了著名的“蒙格馬利”快速冪模演算法：
unsigned Montgomery(unsigned n, unsigned p, unsigned m)
{ // 快速計算 (n ^ e) % m 的值，與power演算法極類似
    unsigned r = n % m; // 這裡的r可不能省
    unsigned k = 1;
    while (p > 1)
    {
        if ((p & 1)!=0)
        {
            k = (k * r) % m; // 直接取模
        }
        r = (r * r) % m; // 同上
        p /= 2;
    }
    return (r * k) % m; // 還是同上
}

上面的程式碼還可以優化。下面是蒙格馬利極速版：

unsigned Montgomery(unsigned n,unsigned p,unsigned m)
{ //快速計算(n^e)%m的值
      unsignedk=1;
      n%=m;
     while(p!=1)
     {
         if(0!=(p&1))k=(k*n)%m;
         n=(n*n)%m;
         p>>=1;
    }
    return(n*k)%m;
}
=====================================================
怎麼判斷一個數是否為素數？

笨蛋的作法： 
bool IsPrime(unsigned n)
{
    if (n<2)
    { //小於2的數即不是合數也不是素數
    throw 0;
    }
    for (unsigned i=2;i<n;++i)
    { //和比它小的所有的數相除，如果都除不盡，證明素數
        if (n%i==0)
        {//除盡了，則是合數
            return false;
        }
    }
    return true;
}

一個數去除以比它的一半還要大的數，一定除不盡，所以還用判斷嗎？？

下面是小學生的做法： 
bool IsPrime(unsigned n)
{
    if (n<2)
    {//小於2的數即不是合數也不是素數
        throw 0;
    }
    for(unsigned i=2;i<n/2+1;++i)
    { // 和比它的一半小數相除，如果都除不盡，證明素數
        if ( 0 == n % i )
        { // 除盡了，合數
            return false;
        }
    }
    return true; // 都沒除盡，素數
}

一個合數必然可以由兩個或多個質數相乘而得到。那麼如果一個數不能被比它的一半小的所有的質數整除，那麼比它一半小的所有的合數也一樣不可能整除它。建立一個素數表是很有用的。

下面是中學生的做法：
bool IsPrime2(unsigned n)
{
    if ( n < 2 )
    { // 小於2的數即不是合數也不是素數
        throw 0;
    }
    static unsigned aPrimeList[] = { // 素數表
        1, 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, 113, 
        193, 241, 257, 337, 353, 401, 433, 449, 577, 593, 641, 
        673, 769, 881, 929, 977, 1009, 1153, 1201, 1217, 1249, 
        1297,1361, 1409, 1489, 1553, 1601, 1697, 1777, 1873, 
        1889, 2017, 2081, 2113, 2129, 2161, 2273, 2417, 2593, 
        2609, 2657, 2689, 2753, 2801, 2833, 2897, 3041, 3089, 
        3121, 3137, 3169, 3217, 3313, 3329, 3361, 3457, 3617, 
        3697, 3761, 3793, 3889, 4001, 4049, 4129, 4177, 4241, 
        4273, 4289, 4337, 4481, 4513, 4561, 4657, 4673, 4721, 
        4801, 4817, 4993, 5009, 5153, 5233, 5281, 5297, 5393, 
        5441, 5521, 5569, 5857, 5953, 6113, 6257, 6337, 6353, 
        6449, 6481, 6529, 6577, 6673, 6689, 6737, 6833, 6961, 
        6977, 7057, 7121, 7297, 7393, 7457, 7489, 7537, 7649, 
        7681, 7793, 7841, 7873, 7937, 8017, 8081, 8161, 8209, 
        8273, 8353, 8369, 8513, 8609, 8641, 8689, 8737, 8753, 
        8849, 8929, 9041, 9137, 9281, 9377, 9473, 9521, 9601, 
        9649, 9697, 9857 
    };
    
    const int nListNum = sizeof(aPrimeList)/sizeof(unsigned);//計算素數表裡元素的個數
    for (unsigned i=2;i<nListNum;++i )
    { 
        if(n/2+1<aPrimeList[i])
        {
            return true;
        }
        if(0==n%aPrimeList[i])
        {
            return false;
        }
    }
    /*由於素數表中元素個數是有限的，那麼對於用素數表判斷不到的數，就只有用笨蛋辦法了*/
    for (unsigned i=aPrimeList[nListNum-1];i<n/2+1;i++ )
    { 
        if (0==n%i)
        { // 除盡了，合數 
            return false;
        }
    }
    return true; 
} 
    還是太糟了，我們現在要做的對於大型素數的判斷，那個素數表倒頂個P用！當然，我們可以利用動態的素數表來進行優化，這就是大學生的做法了。但是動態生成素數表的策略又複雜又沒有效率，所以我們還是直接跳躍到專家的做法吧：
    根據上面講到的費馬小定理，對於兩個互質的素數N和P，必有：N^(P-1)%P=1 
    那麼我們通過這個性質來判斷素數吧，當然，你會擔心當P很大的時候乘方會很麻煩。不用擔心！我們上面不是有個快速的冪模演算法麼？好好的利用蒙格馬利這位大數學家為我們帶來的快樂吧！

演算法思路是這樣的： 
    對於N，從素數表中取出任意的素數對其進行費馬測試，如果取了很多個素數，N仍未測試失敗，那麼則認為N是素數。當然，測試次數越多越準確，但一般來講 50次就足夠了。另外，預先用“小學生”的演算法構造一個包括500個素數的陣列，先對Q進行整除測試，將會大大提高通過率，方法如下:
bool IsPrime3(unsigned n)
{
    if ( n < 2 )
    { // 小於2的數即不是合數也不是素數
        throw 0;
    }
    static unsigned aPrimeList[] = {
        2, 3, 5, 7, 11, 17, 19, 23, 29, 31, 41,
        43, 47, 53, 59, 67, 71, 73, 79, 83, 89, 97
    };
    const int nListNum = sizeof(aPrimeList) / sizeof(unsigned);
    for (int i=0;i<nListNum;++i)
    { // 按照素數表中的數對當前素數進行判斷
        if (1!=Montgomery(aPrimeList[i],n-1,n)) // 蒙格馬利演算法
        {
            return false;
        }
    }
    return true;
}

    OK，這就專家的作法了。 
    等等，什麼？好像有點怪，看一下這個數29341，它等於13 * 37 * 61，顯然是一個合數，但是竟通過了測試！！哦，抱歉，我忘了在素數表中加入13，37，61這三個數，我其實是故意的，我只是想說明並費馬測試並不完全可靠。
    現在我們發現了重要的一點，費馬定理是素數的必要條件而非充分條件。這種不是素數，但又能通過費馬測試的數字還有不少，數學上把它們稱為卡爾麥克數，現在數學家們已經找到所有10 ^ 16以內的卡爾麥克數，最大的一個是9585921133193329。我們必須尋找更為有效的測試方法。數學家們通過對費馬小定理的研究，並加以擴充套件，總結出了多種快速有效的素數測試方法，目前最快的演算法是拉賓米勒測試演算法，下面介紹拉賓米勒測試。
================================================================
拉賓米勒測試

    拉賓米勒測試是一個不確定的演算法，只能從概率意義上判定一個數可能是素數，但並不能確保。演算法 流程如下:
    1.選擇T個隨機數A，並且有A<N成立。
    2.找到R和M，使得N=2*R*M+1成立。
    快速得到R和M的方式：N用二進位制數B來表示，令C=B-1。因為N為奇數（素數都是奇數），所以C的最低位為0，從C的最低位的0開始向高位統計，一直到遇到第一個1。這時0的個數即為R，M為B右移R位的值 。
    3.如果A^M%N=1，則通過A對於N的測試，然後進行下一個A的測試
    4.如果A^M%N!=1，那麼令i由0迭代至R，進行下面的測試
    5.如果A^((2^i)*M)%N=N-1則通過A對於N的測試，否則進行下一個i的測試 
    6.如果i=r，且尚未通過測試，則此A對於N的測試失敗，說明N為合數。
    7.進行下一個A對N的測試，直到測試完指定個數的A

    通過驗證得知，當T為素數，並且A是平均分佈的隨機數，那麼測試有效率為1 / ( 4 ^ T )。如果T > 8那麼測試失誤的機率就會小於10^(-5)，這對於一般的應用是足夠了。如果需要求的素數極大，或著要求更高的保障度，可以適當調高T的值。下面是程式碼：

bool RabbinMillerTest( unsigned n ) 
{
    if (n<2)
    { // 小於2的數即不是合數也不是素數
        throw 0;
    }

    const unsigned nPrimeListSize=sizeof(g_aPrimeList)/sizeof(unsigned);//求素數表元素個數
    for(int i=0;i<nPrimeListSize;++i)
    {// 按照素數表中的數對當前素數進行判斷
        if (n/2+1<=g_aPrimeList[i])
        {// 如果已經小於當前素數表的數，則一定是素數
            return true;
        }
        if (0==n%g_aPrimeList[i])
        {// 餘數為0則說明一定不是素數
            return false;
        }
    }
    // 找到r和m，使得n = 2^r * m + 1;
    int r = 0, m = n - 1; // ( n - 1 ) 一定是合數
    while ( 0 == ( m & 1 ) )
    {
        m >>= 1; // 右移一位
        r++; // 統計右移的次數
    }
    const unsigned nTestCnt = 8; // 表示進行測試的次數
    for ( unsigned i = 0; i < nTestCnt; ++i )
    { // 利用隨機數進行測試，
        int a = g_aPrimeList[ rand() % nPrimeListSize ];
        if ( 1 != Montgomery( a, m, n ) )
        {
            int j = 0;
            int e = 1;
            for ( ; j < r; ++j )
            {
                if ( n - 1 == Montgomery( a, m * e, n ) ) 
                {
                    break;
                }
                e <<= 1;
            }
            if (j == r)
            {
                return false;
            }
        }
    }
    return true;
}