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

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

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

3、把N提出來一個，NP就成了你N*(N(P-1)，那么(NP-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的公倍數(shù)，小學(xué)知識了_又因為前提是N與P互質(zhì)，而互質(zhì)數(shù)的最小公倍數(shù)為它們的乘積，所以一定存在正整數(shù)M使得等式成立：N*(N(P-1)-1)=M*N*P兩邊約去N，化簡之：N(P-1)-1=M*P因為M是整數(shù)，顯然：(N(P-1)-1)%P=0即：N(P-1)%P=1=積模分解公式先有一個引理，如果有：X%Z=0，即X能被Z整除，則有：

4、(X+Y)%Z=Y%Z這個不用證了吧.設(shè)有X、Y和Z三個正整數(shù)，則必有：(X*Y)%Z=(X%Z)*(Y%Z)%Z想了很長時間才證出來，要分情況討論才行：1.當(dāng)X和Y都比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ù)倍暈，又來了。概據(jù)引理，（3）式可化簡為：(A*B)%Z又因為：A=X%Z，B=Y%Z，代入上面的

5、式子，就成了原式了。2.當(dāng)X比Z大而Y比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，代入上式，即得到原式。同理，當(dāng)X比Z小而Y比Z大時，原式也成立。3.當(dāng)X比Z小，且Y也比Z小時，X=X%Z，Y=Y%Z，所以原式成立。=快速計算乘方的算法如計算213，則傳統(tǒng)做法需要進(jìn)行12次乘法。/*計算np*/unsigned power(unsigned n,unsigned p) for(int i=0;i<p;i+) n*=n; return n;該死的乘法，是時候優(yōu)化一下了！把2*2的結(jié)果保

6、存起來看看，是不是成了：4*4*4*4*4*4*2 再把4*4的結(jié)果保存起來：16*16*16*2 一共5次運算，分別是2*2、4*4和16*16*16*2這樣分析，我們算法因該是只需要計算一半都不到的乘法了。為了講清這個算法，再舉一個例子27：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 實際上最后一個括號里的2 * 2是這回又剩下的，那么，稍后再單獨乘上它 現(xiàn)在指數(shù)已經(jīng)為1了，可以計算最終結(jié)果了：16*4*2=128

7、優(yōu)化后的算法如下：unsigned Power(unsigned n,unsigned p) unsigned main=n; /用main保存結(jié)果unsigned odd=1; /odd用來計算“剩下的”乘積while (p>1) /一直計算，直到指數(shù)小于或等于1 if(p%2)!=0) / 如果指數(shù)p是奇數(shù)，則說明計算后會剩一個多余的數(shù)，那么在這里把它乘到結(jié)果中 odd*=main; /把“剩下的”乘起來 main*=main; /主體乘方 p/=2; /指數(shù)除以2 return main*odd; /最后把主體和“剩下的”乘起來作為結(jié)果夠完美了嗎？不，還不夠！看出來了嗎？main是

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

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

10、將這個過程優(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 = (r * r) % m; / 同上 p /= 2; return (r * k) % m; / 還是同上上面的代碼還可以優(yōu)化。下面是蒙格馬

11、利極速版：unsigned Montgomery(unsigned n,unsigned p,unsigned m) /快速計算(ne)%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ù)也不是素數(shù) throw 0; for (unsigned i=2;i<n;+i) /和比它小的所有的數(shù)相除，如果都除不盡

12、，證明素數(shù) if (n%i=0) /除盡了，則是合數(shù) return false; return true;一個數(shù)去除以比它的一半還要大的數(shù)，一定除不盡，所以還用判斷嗎？下面是小學(xué)生的做法： 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 false; return true; / 都沒除盡，素數(shù)一個合數(shù)必然可以由兩個或多個質(zhì)數(shù)相乘而得到。那

13、么如果一個數(shù)不能被比它的一半小的所有的質(zhì)數(shù)整除，那么比它一半小的所有的合數(shù)也一樣不可能整除它。建立一個素數(shù)表是很有用的。下面是中學(xué)生的做法：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, 89, 97, 113, 193, 241, 257, 337, 35

14、3, 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, 361

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

16、3, 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<

17、;nListNum;+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用！當(dāng)然，我們可以利用動態(tài)的素數(shù)表來進(jìn)行優(yōu)化，這就是大學(xué)生的做法了。但是動

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

19、，方法如下：bool IsPrime3(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ù)對當(dāng)前素數(shù)進(jìn)行判斷

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

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

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