費(fèi)馬小定理素?cái)?shù)判定蒙哥馬利算法

1、費(fèi)馬小定理素?cái)?shù)判定蒙哥馬利算法(強(qiáng)烈推薦)2009-11-07 12:42費(fèi)馬小定理素?cái)?shù)判定蒙哥馬利算法約定：x%y為x取模y,即x除以y所得的余數(shù)，當(dāng)x<y時(shí)，x%y=x ,所有取模的運(yùn)算對(duì)象都為整數(shù)。xAy 表示x的y次方。乘方運(yùn)算的優(yōu)先級(jí)高于乘除和取模，加減的優(yōu)先級(jí)最低。見到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的一樣，靠！復(fù)習(xí)一下小學(xué)數(shù)學(xué)公因數(shù)：兩個(gè)不同的自然數(shù)A和B,若有自然數(shù)C可以整除A也可以整除B,那么C就是A和B的公因數(shù)。

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

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

4、模分解公式先有一個(gè)引理，如果有：X%Z=0 ,即X能被Z整除，則有：(X+Y)%Z=Y%Z這個(gè)不用證了吧設(shè)有X、丫和Z三個(gè)正整數(shù)，則必有：(X*Y)%Z=(X%Z)*(Y%Z)%Z想了很長(zhǎng)時(shí)間才證出來，要分情況討論才行：1.當(dāng)X和丫都比Z大時(shí)，必有整數(shù) A和B使下面的等式成立：X=Z*I+A ( 1)Y=Z*J+B (2)不用多說了吧，這是除模運(yùn)算的性質(zhì)！將(1)和(2)代入(X*Y)modZ 得：(Z*I+A)(Z* J+B)%Z乘開，再把前三項(xiàng)的 Z提一個(gè)出來，變形為：(Z*(Z*I* J+I*A+I*B)+A*B)%Z(3)因?yàn)閆*(Z*I* J+I*A+I*B)是Z的整數(shù)倍暈，又來了。

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

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

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

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

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

10、可以把丫個(gè)X相乘再取模的過程分解開來，比如： (17人25)%29 則可分解為：(17 *17) % 29 * ( 17*17) % 29 *如果用上面的代碼將這個(gè)過程優(yōu)化，那么我們就得到了著名的“蒙格馬利”快速募模算法：unsigned Montgomery(unsigned n, unsigned p, unsigned m) /快速計(jì)算(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) /快速計(jì)算(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;怎么判斷一個(gè)數(shù)是否為素?cái)?shù)？笨蛋的作法：bool IsPrime(unsigned n)if (n<2) /小于2的數(shù)即不是合數(shù)也不

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

13、e;return true; /都沒除盡，素?cái)?shù)一個(gè)合數(shù)必然可以由兩個(gè)或多個(gè)質(zhì)數(shù)相乘而得到。那么如果一個(gè)數(shù)不能被比它的一半小的所有的質(zhì)數(shù)整除，那么比它一半小的所有的合數(shù) 也一樣不可能整除它。建立一個(gè)素?cái)?shù)表是很有用的。下面是中學(xué)生的做法：bool IsPrime2(unsigned n)if ( n < 2 ) 小于2的數(shù)即不是合數(shù)也不是素?cái)?shù)throw 0;static unsigned aPrimeList = / 素?cái)?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);計(jì)算素?cái)?shù)表里元素的個(gè)數(shù)for (unsigned i=2;i<nLi

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

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

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

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

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

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