1.約定 x%y為x取模y，即x除以y所得的余數，當x<y時，x%y=x，所有取模的運算對象都為整數。 x^y表示x的y次方。乘方運算的優先級高於乘除和取模，加減的優先級最低。 見到x^y/z這樣，就先算乘方，再算除法。www.2cto.com A/B，稱為A除以B，也稱為B除A。 若A%B=0，即稱為A可以被B整除，也稱B可以整除A。 A*B表示A乘以B或稱A乘B，B乘A，B乘以A……都一樣。 復習一下小學數學 公因數：兩個不同的自然數A和B，若有自然數C可以整除A也可以整除B，那麼C就是A和B的公因數。 公倍數：兩個不同的自然數A和B，若有自然數C可以被A整除也可以被B整除，那麼C就是A和B的公倍數。 互質數：兩個不同的自然數，它們只有一個公因數1，則稱它們互質。 費馬是法國數學家，又譯“費爾馬”，此人巨牛，他的簡介請看下面。不看不知道，一看嚇一跳。 2.費馬小定理： 有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 ============================================ 3.積模分解公式 先有一個引理，如果有：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，所以原式成立。 ===================================================== 4.快速計算乘方的算法 如計算2^13，則傳統做法需要進行12次乘法。 [cpp] /*計算n^p*/ unsigned power(unsigned n,unsigned p) { for(int i=0;i<p;i++) n*=n; return n; } /*計算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 優化後的算法如下： [cpp] 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; //最後把主體和“剩下的”乘起來作為結果 } 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! 完美版： [cpp] 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; // 最後把主體和“剩下的”乘起來作為結果 } 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; // 最後把主體和“剩下的”乘起來作為結果 } ======================================================== 5."蒙格馬利”快速冪模算法 後面我們會用到這樣一種運算：(X^Y)%Z。但問題是當X和Y很大時，只有32位的整型變量如何能夠有效的計算出結果？ 考慮上面那份最終的優化代碼和再上面提到過的積模分解公式，我想你也許會猛拍一下腦門，吸口氣說：“哦，我懂了！”。 下面的講解是給尚沒有做出這樣動作的同學們准備的： X^Y可以看作Y個X相乘，即然有積模分解公式，那麼我們就可以把Y個X相乘再取模的過程分解開來，比如：(17^25)%29則可分解為：( ( 17 * 17 ) % 29 * ( 17 * 17 ) % 29 * …… 如果用上面的代碼將這個過程優化，那麼我們就得到了著名的“蒙格馬利”快速冪模算法： [cpp] 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 的值，與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; // 還是同上 } 上面的代碼還可以優化。下面是蒙格馬利極速版： [cpp] unsigned Montgomery(unsigned n,unsigned p,unsigned m) { //快速計算(n^p)%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; } unsigned Montgomery(unsigned n,unsigned p,unsigned m) { //快速計算(n^p)%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; } ===================================================== 6.怎麼判斷一個數是否為素數？ 1)笨蛋的作法： [cpp] 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;++i) { //和比它小的所有的數相除，如果都除不盡，證明素數 if (n%i==0) { //除盡了，則是合數 return false; } } return true; }一個數去除以比它的一半還要大的數，一定除不盡，所以還用判斷嗎？？ 2)下面是小學生的做法： [cpp] 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 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; // 都沒除盡，素數 } 一個合數必然可以由兩個或多個質數相乘而得到。那麼如果一個數不能被比它的一半小的所有的質數整除，那麼比它一半小的所有的合數也一樣不可能整除它。建立一個素數表是很有用的。 3)下面是中學生的做法： [cpp] 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; } 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進行整除測試，將會大大提高通過率，方法如下: 6)下面是專家的做法： [cpp] 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; } 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。我們必須尋找更為有效的測試方法。數學家們通過對費馬小定理的研究，並加以擴展，總結出了多種快速有效的素數測試方法，目前最快的算法是拉賓米勒測試算法，下面介紹拉賓米勒測試。 ================================================================ 7.拉賓米勒測試 拉賓米勒測試是一個不確定的算法，只能從概率意義上判定一個數可能是素數，但並不能確保。算法流程如下: 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的值。 下面是代碼： [cpp] 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; } 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; }
