費馬小定理素數(shù)判定蒙哥馬利算法

1、費馬小定理素數(shù)判定蒙哥馬利算法(強烈推薦)2009-11-07 12:42費馬小定理素數(shù)判定蒙哥馬利算法約定：x%y為x取模y,即x除以y所得的余數(shù)，當x<y時，x%y=x ,所有取模的運算對象都為整數(shù)。xAy 表示x的y次方。乘方運算的優(yōu)先級高于乘除和取模，加減的優(yōu)先級最低。見到xAy/z這樣，就先算乘方，再算除法。A/B ,稱為A除以B,也稱為B除A。若A%B=0 ,即稱為A可以被B整除，也稱B可以整除AoA*B表示A乘以B或稱A乘B, B乘A, B乘以A都TMD的一樣，靠！復習一下小學數(shù)學公因數(shù)：兩個不同的自然數(shù)A和B,若有自然數(shù)C可以整除A也可以整除B,那么C就是A和B的公因數(shù)。

2、公倍數(shù)：兩個不同的自然數(shù)A和B,若有自然數(shù)C可以被A整除也可以被B整除，那么C就是A和B的公倍數(shù)?；ベ|(zhì)數(shù)：兩個不同的自然數(shù)，它們只有一個公因數(shù)1,則稱它們互質(zhì)。費馬是法國數(shù)學家，又譯“費爾馬”，此人巨牛，他的簡介請看下面。不看不知道，一看嚇一跳。費馬小定理：有N為任意正整數(shù)，P為素數(shù)，且N不能被P整除(顯然N和P互質(zhì))，則有：NAP%P=N( 即：N的P次方除以 P的余數(shù)是 N)但是我查了很多資料見到的公式都是這個樣子：(NA(P-1)%P=1后來分析了一下，兩個式子其實是一樣的，可以互相變形得到，原式可化為：(NAP-N)%P=0( 即：N的P次方減N可以被P整除，因為由費馬小定理知道N的P

3、次方除以P的余數(shù)是N)把 N 提出來一個，NAP 就成了你 N*(NA(P-1)，那么(NAP-N)%P=0可化為：(N*(NA(P-1)-1)%P=0請注意上式，含義是：N*(NA(P-1)-1) 可以被P整除又因為N*(NA(P-1)-1) 必能整除N (這不費話么！)所以，N*(NA(P-1)-1) 是N和P的公倍數(shù)，小學知識了 a_a又因為前提是N與P互質(zhì)，而互質(zhì)數(shù)的最小公倍數(shù)為它們的乘積，所以一定存在正整數(shù)M使得等式成立：N*(NA(P-1)-1)=M*N*P兩邊約去N ,化簡之：NA(P-1)-1=m*p因為M是整數(shù)，顯然：(NA(P-1)-1)%P=0即：NA(P-1)%P=1積

4、模分解公式先有一個引理，如果有：X%Z=0 ,即X能被Z整除，則有：(X+Y)%Z=Y%Z這個不用證了吧設(shè)有X、丫和Z三個正整數(shù)，則必有：(X*Y)%Z=(X%Z)*(Y%Z)%Z想了很長時間才證出來，要分情況討論才行：1.當X和丫都比Z大時，必有整數(shù) A和B使下面的等式成立：X=Z*I+A ( 1)Y=Z*J+B (2)不用多說了吧，這是除模運算的性質(zhì)！將(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的整數(shù)倍暈，又來了。

5、概據(jù)引理，(3)式可化簡為：(A*B)%Z又因為：A=X%Z , B=Y%Z ,代入上面的式子，就成了原式了。2 .當X比Z大而丫比Z小時，一樣的轉(zhuǎn)化：X=Z*I+A代入(X*Y)%Z得：(Z*I*Y+A *Y)%Z根據(jù)引理，轉(zhuǎn)化得：(A*Y)%Z因為A=X%Z ,又因為Y=Y%Z ,代入上式，即得到原式。同理，當X比Z小而丫比Z大時，原式也成立。3 .當X比Z小，且 丫也比Z小時，X=X%Z , Y=Y%Z ,所以原式成立?？焖儆嬎愠朔降乃惴ㄈ缬嬎?A13 ,則傳統(tǒng)做法需要進行12次乘法。/* 計算 nAp*/unsigned power(unsigned "unsigned p)

6、for(int i=0;i<p;i+) n*=n;return n;該死的乘法，是時候優(yōu)化一下了！把 2*2的結(jié)果保存起來看看，是不是成了：4*4*4*4*4*4*2再把4*4的結(jié)果保存起來：16*16*16*2一共5次運算，分別是2*2、4*4和16*16*16*2這樣分析，我們算法因該是只需要計算一半都不到的乘法了。為了講清這個算法，再舉一個例子2A7 : 2*2*2*2*2*2*2兩兩分開:(2*2)*(2*2)*(2*2)*2如果用2*2來計算，那么指數(shù)就可以除以2了，不過剩了一個，稍后再單獨乘上它。再次兩兩分開，指數(shù)除以 2: (2*2)*(2*2)*(2*2)*2實際上最后一

7、個括號里的 2*2是這回又剩下的，那么，稍后再單獨乘上它現(xiàn)在指數(shù)已經(jīng)為1 了，可以計算最終結(jié)果了：16*4*2=128優(yōu)化后的算法如下：unsigned Power(unsigned "unsigned p)unsigned main=n; 用 main 保存結(jié)果unsigned odd=1;/odd 用來計算“剩下的"乘積while (p>1)/ 一直計算，直到指數(shù)小于或等于1if(p%2)!=0)/如果指數(shù)p是奇數(shù)，則說明計算后會剩一個多余的數(shù)，那么在這里把它乘到結(jié)果中 odd*=main;/ 把"剩下的"乘起來main*=main; /主體乘

8、方p/=2; /指數(shù)除以2return main* odd; /最后把主體和“剩下的"乘起來作為結(jié)果夠完美了嗎？不，還不夠！看出來了嗎？ main是沒有必要的，并且我們可以有更快的代碼來判斷奇數(shù)。要知道除法或取模運算的效率很低，所以我們可以利用偶數(shù)的一個性質(zhì)來優(yōu)化代碼，那就是偶數(shù)的二進制表示法中的最低位一定為0!完美版：unsigned Power(unsigned n, unsigned p) /計算n的p次方unsigned odd = 1; /oddk 用來計算“剩下的"乘積while (p > 1) / 一直計算到指數(shù)小于或等于 1if ( p & 1

9、 )!=0) /判斷p是否奇數(shù)，偶數(shù)的最低位必為0odd *= n; / 若odd為奇數(shù)，則把“剩下的"乘起來n *= n; /主體乘方p /= 2; /指數(shù)除以2return n * odd; /最后把主體和“剩下的"乘起來作為結(jié)果蒙格馬利”快速募模算法后面我們會用到這樣一種運算：(XAY)%Z問題是當X和丫很大時，只有32位的整型變量如何能夠有效的計算出結(jié)果？考慮上面那份最終的優(yōu)化代碼和再上面提到過的積模分解公式，我想你也許會猛拍一下腦門，吸口氣說：“哦，我懂了！”。下面的講解是給尚沒有做出這樣動作的同學們準備的。XAY可以看作丫個X相乘，即然有積模分解公式，那么我們就

10、可以把丫個X相乘再取模的過程分解開來，比如： (17人25)%29 則可分解為：(17 *17) % 29 * ( 17*17) % 29 *如果用上面的代碼將這個過程優(yōu)化，那么我們就得到了著名的“蒙格馬利”快速募模算法：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

11、 = (r * r) % m; / 同上p /= 2;return (r * k) % m; / 還是同上上面的代碼還可以優(yōu)化。下面是蒙格馬利極速版：unsigned Montgomery(unsigned n,unsigned p,unsigned m) /快速計算(Me)%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;怎么判斷一個數(shù)是否為素數(shù)？笨蛋的作法：bool IsPrime(unsigned n)if (n<2) /小于2的數(shù)即不是合數(shù)也不

12、是素數(shù)throw 0;for (unsigned i=2;i<n;+i) /和比它小的所有的數(shù)相除，如果都除不盡，證明素數(shù)if (n%i=0)/除盡了,則是合數(shù)return false;return true;一個數(shù)去除以比它的一半還要大的數(shù)，一定除不盡，所以還用判斷嗎？ ？下面是小學生的做法：bool IsPrime(unsigned n)if (n<2)/小于2的數(shù)即不是合數(shù)也不是素數(shù)throw 0;for(unsigned i=2;i<n/2+1;+i) /和比它的一半小數(shù)相除，如果都除不盡，證明素數(shù)if ( 0 = n % i ) /除盡了，合數(shù)return fals

13、e;return true; /都沒除盡，素數(shù)一個合數(shù)必然可以由兩個或多個質(zhì)數(shù)相乘而得到。那么如果一個數(shù)不能被比它的一半小的所有的質(zhì)數(shù)整除，那么比它一半小的所有的合數(shù) 也一樣不可能整除它。建立一個素數(shù)表是很有用的。下面是中學生的做法：bool IsPrime2(unsigned n)if ( n < 2 ) 小于2的數(shù)即不是合數(shù)也不是素數(shù)throw 0;static unsigned aPrimeList = / 素數(shù)表1,2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31,37, 41,43, 47, 53, 59, 61,67, 71, 73, 79, 83

14、, 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, 33

15、13, 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,

16、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);計算素數(shù)表里元素的個數(shù)for (unsigned i=2;i<nLi

17、stNum;+i )if(n/2+1<aPrimeListi)return true;if(0=n%aPrimeListi)return false;/*由于素數(shù)表中元素個數(shù)是有限的，那么對于用素數(shù)表判斷不到的數(shù)，就只有用笨蛋辦法了*/for (unsigned i=aPrimeListnListNum-1;i<n/2+1;i+)if (0=n%i) 除盡了，合數(shù) return false;return true;還是太糟了，我們現(xiàn)在要做的對于大型素數(shù)的判斷，那個素數(shù)表倒頂個P用！當然，我們可以利用動態(tài)的素數(shù)表來進行優(yōu)化，這就是大學生的做法了。但是動態(tài)生成素數(shù)表的策略又復雜又沒有效

18、率，所以我們還是直接跳躍到專家的做法吧：根據(jù)上面講到的費馬小定理，對于兩個互質(zhì)的素數(shù) N和P,必有：NA(P-1)%P=1那么我們通過這個性質(zhì)來判斷素數(shù)吧，當然，你會擔心當P很大的時候乘方會很麻煩。不用擔心！我們上面不是有個快速的嘉模算法么？好好的利用蒙格馬利這位大數(shù)學家為我們帶來的快樂吧！ 算法思路是這樣的：對于N ,從素數(shù)表中取出任意的素數(shù)對其進行費馬測試，如果取了很多個素數(shù)，N仍未測試失敗，那么則認為 N是素數(shù)。當然，測試次數(shù)越多越準確，但一般來講 50次就足夠了。另外，預先用“小學生”的算法構(gòu)造一個包括500個素數(shù)的數(shù)組，先對 Q進行整除測試，將會大大提高通過率，方法如下：bool I

19、sPrime3(unsigned n)if ( n < 2 ) /小于2的數(shù)即不是合數(shù)也不是素數(shù) 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) /按照素數(shù)表中的數(shù)對當前素數(shù)進行判斷if (1!=Montgomery(aPri

20、meListi,n-1,n) /蒙格馬利算法return false;return true;OK,這就專家的作法了。等等，什么？好像有點怪，看一下這個數(shù)29341 ,它等于13 * 37 * 61 ,顯然是一個合數(shù)，但是竟通過了測試！哦，抱歉，我忘了在素數(shù)表中加入13, 37, 61這三個數(shù)，我其實是故意的，我只是想說明并費馬測試并不完全可靠?，F(xiàn)在我們發(fā)現(xiàn)了重要的一點，費馬定理是素數(shù)的必要條件而非充分條件。這種不是素數(shù)，但又能通過費馬測試的數(shù)字還有不少，數(shù)學上把它們稱為卡爾麥克數(shù)，現(xiàn)在數(shù)學家們已經(jīng)找到所有10人16以內(nèi)的卡爾麥克數(shù)，最大的一個是 9585921133193329。我們必須尋找

21、更為有效的測試方法。數(shù)學家們通過對費馬小定理的研究，并加以擴展，總結(jié)出了多種快速有效的素數(shù)測試方法，目前最快的算法是拉賓米 勒測試算法，下面介紹拉賓米勒測試。拉賓米勒測試拉賓米勒測試是一個不確定的算法，只能從概率意義上判定一個數(shù)可能是素數(shù)，但并不能確保。算法流程如下1 .選擇T個隨機數(shù)A,并且有A<N成立。2 .找到R和M ,使得N=2*R*M+1 成立?？焖俚玫絉和M的方式：N用二進制數(shù)B來表示，令C=B-1 o因為N為奇數(shù)(素數(shù)都是奇數(shù))，所以C的最低位為0,從C的最低 位的0開始向高位統(tǒng)計，一直到遇到第一個1。這時0的個數(shù)即為R, M為B右移R位的值。3 .如果AAM%N=1,則通

22、過A對于N的測試，然后進行下一個A的測試4 .如果AAM%N!=1,那么令i由0迭代至R,進行下面的測試5 .如果AA(2Ai)*M)%N=N-1 則通過A對于N的測試，否則進行下一個i的測試6 .如果i=r ,且尚未通過測試，則此A對于N的測試失敗，說明 N為合數(shù)。7 .進行下一個A對N的測試，直到測試完指定個數(shù)的A通過驗證得知，當 T為素數(shù)，并且A是平均分布的隨機數(shù)，那么測試有效率為1 / ( 4 a T )。如果T > 8那么測試失誤的機率就會小于10A(-5),這對于一般的應用是足夠了。如果需要求的素數(shù)極大，或著要求更高的保障度，可以適當調(diào)高T的值。下面是代碼：bool RabbinMillerTest