1. 程式人生 > >米勒-拉賓(MillerRabbin)素性測試演算法詳解

米勒-拉賓(MillerRabbin)素性測試演算法詳解

寫在前面

網上有很多關於米勒拉賓素性測試演算法的部落格  但是大多數都是轉載,或者只有模板程式碼沒有分析講解的,甚至還有的分析的都是錯的。花了一早上,借鑑了幾十篇部落格,總算是把這個演算法理解了差不多,並且詳細整理了一下我的理解。

講解很細,篇幅較長,要是想看,準備好耐心。

個人理解,如有不對,歡迎指正。

 

米勒-拉賓(MillerRabbin)素性測試演算法

 

米勒-拉賓(MillerRabbin)素性測試演算法是一個高效判斷素數的方法

 

首先 在學習米勒拉賓素性測試涉及到兩個定理  先放出來  下文會詳細解釋

1、費馬小定理: 如果p為質數           a^{p-1}\equiv 1%p (在mod p 的情況下)

2、對於任意一個小於p的正整數x,發現1(模p)的非平凡平方根存在,則說明p是合數。

首先說  費馬小定理

費馬素性測試

關於 費馬小定理   他的證明 這裡就不贅述了  感興趣可以自行百度

我們要用的是費馬小定理的逆推

費馬小定理說, 如果p為質數    則滿足 a^{p-1}\equiv 1%p (在mod p 的情況下)

那麼是不是意味著  如果任意一個數p   滿足 a^{p-1}\equiv 1%p

(mod p )  就可以斷定p是質數了呢?

如果這樣來檢測素數   這個方法就叫做費馬素性測試      但是答案是  不可以這樣!

  要注意  p是質數  一定滿足費馬小定理   但是  滿足費馬小定理的數  並不能證明就是質數

有些數字 滿足費馬小定理,但是並不是質數  他們叫做偽質數 Carmichael(偽質數的個數是無窮的  文末附有偽質數表哦)

所以說 如果我們用費馬檢測去判斷素數是會出錯的

那麼這個出錯的概率是多大呢    在一定範圍內 有多少個數是偽素數呢?  這個和a的選取有關係

當a = 2時,前十億個正整數中能滿足費馬小定理且不是素數的只有5597個  ,也就是說 出錯的概率是0.011%

雖然這個錯誤率不高 但是還是存在那麼多的偽素數  費馬檢測無法檢測出來

我需要一個新的方法 來判斷這些偽素數

這個時候也就是米勒拉賓演算法的核心部分了  二次探測

二次探測

這裡就要用到一開始放出來的那個定理二

2、對於任意一個小於p的正整數x,發現1(模p)的非平凡平方根存在,則說明p是合數。

乍一看是不是根本看不懂這個定理說的啥意思

其實翻譯下來就是下面這樣

如果p是一個素數,0 < x< p,   則方程 x^{2} ≡ 1(mod p) 的解為 x=1 ,x= p-1

反之 如果  x^ 2 ≡ 1(mod p)  的解不是 x=1 ,x= p-1   那 p 就不是素數

(我們把 1和-1稱為 平凡根   其他的根都稱為非平凡根

而 如果一個數x 的平方 模 p =1 或-1   我們就稱這個數是1 或 -1 模p 的非平凡平方根

然而 x^{2} ≡ -1 (mod p)  那x就是虛數了 不在我們討論範圍內  所以  我們只用關心 x^{2} ≡ 1(mod p)  

也就是隻用關心 1(模p)的非平凡平方根  而如果p是素數 那麼 他的根只可能是 1或者 p-1

若根不是1或者p-1 那他就是個合數     關於證明 在下邊)

我們用這條性質 進行二次探測

關於這條性質的證明: (如果不想看證明直接跳過 往下看

--------------------------------證明-----------------------------------------

x^{2} ≡ 1(mod p)   可以移項  變成    x^{2}-1 ≡ 0 (mod p)

那麼就是  ( x - 1 )( x + 1 ) ≡ 0(mod p)     就是說  ( x - 1 )( x + 1 ) % p = 0

那如果p是個質數

因為p是質數     所以 p不可能  %( x - 1 )=0   也不可能 %( x + 1 )=0

既然%他倆都不可能為0  也就是他倆都不是p的因子(x<p)   那他倆的乘積也不可能%p = 0

x^{2} ≡ 1(mod p)  這怎麼成立呢

只有一個答案  ( x - 1 )( x + 1 )本身就=0   那麼也就是  x=1 或 p-1

由此可證明  p為素數的時候 x^{2} ≡ 1(mod p)   只有兩個解 x=1 或 p-1

--------------------------------證畢--------------------------------------------

有了這條性質 

那麼    滿足x^ 2 ≡ 1(mod p)  的 x 值   只可能等於  1 或 p-1     如果不是這倆   那就不是素數 (x<p)

那麼我們去怎樣驗證呢  總不能把比 p 小的數一個一個驗證吧  那樣也太費時間了吧

所以 米勒拉賓演算法就提供了一個 快速的辦法(但不準確~~後面會講)

我們轉換個角度看問題

既然   我們遍歷所有 x 讓它 mod  p  結果有很多種   我們需要的只是那些結果為1的 x

那我們反過來   把所有結果為 1 的 x 拿出來       看他們 是不是等於   p - 1 或 1   

 

寫這裡 你可能要說了  what?  怎麼根據結果找這些 x    我不經過計算怎麼知道結果是不是1   ?

彆著急   可以的  首先要明白   我們要找的 這些 x     都是滿足 x ^ 2 ≡ 1 (mod p) (別忘記了~~)

那麼 再看下剛才的費馬素性檢測

能通過費馬檢測的數    一定滿足     a^{p-1}\equiv 1%p (mod  p)

而 p - 1是一個偶數

(如果p是個偶數   我們不需要任何複雜的演算法  因為偶數只有2是質數  直接判斷即可  所以只有p是奇數 我們才進來判斷   p是奇數  p - 1  自然是偶數了)

然後 來看個很好理解的性質

任何一個偶數 都可以寫成   2^{r}*d的形式      比如 6 可以寫成  2^1 * 3

所以  p-1 也可以寫成 2^{r}*d   

那麼  a^{p-1} 就等於 a^{2^{r}*d}      也就可以寫成 (a^{d})^{2^{r}}

那麼式子就變成了  (a^{d})^{2^{r}}\equiv 1(mod p)(mod  p)

我們把  a^{d}  看作 k       那麼式子就是      (k)^{2^{r}} \equiv 1(mod  p) 

 

為什麼把費馬小定理化簡成這個樣子呢?  我們來看看我們要找的 x    x是滿足 x ^ 2 ≡ 1 (mod p) 解為 x=1 ,x= p-1的

而  (k)^{2^{r}}的計算過程      是不是 (k*k)^{2^{r}-1}   ....(k*k*k)^{2^{r}-2}      (k*k*k*k)^{2^{r}-3}................

2^{r}的變化過程中  一定會存在    (k*k*...*k)^{2}    那麼  此時  這些k相乘  就是我們要找的 x

這就是滿足   x ^ 2 ≡ 1 (mod p) 的那些 x        這個x是存在於 (k)^{2^{r}} 的變化過程中的

所以我們只需要在2^{r}的變化過程中不停的判斷  只要它等於1 或者 p - 1  就說明他是素數!

當然 這裡還有一個重點    如果 x  在這個過程中一直不等於  1或 p-1    就能說明他不是素數嗎

不能!   思路不要暈   我們要徹底明白米勒拉賓的核心  

核心是 根據二次探測定理   我們要找滿足 x ^ 2 ≡ 1 (mod p) 的所有  x 

要看這些x是不是全部為 1 或 p-1   

但是挨個遍歷太慢了

然後我們發現  滿足x ^ 2 ≡ 1 (mod p)  的x    都存在於  (k)^{2^{r}} 裡面

然後我們去 (k)^{2^{r}} 裡面找 x    但是 (k)^{2^{r}}是和 k 有關的

嚴格來說    你想判斷一個數不是素數   需要找遍所有滿足條件的 x 

而想找遍 x    你必須找遍所有k的取值

明白了這個 你可能會說 找遍所有k 肯定超時啊   扯了半天 這不還是判斷不了嗎。。。

別急 確實 米勒拉賓不能百分百的判斷出一個數是不是素數

但是  我們只嘗試幾個 k 值    去找x    如果x取值都不是 1或 p-1  雖然理論上不能斷定他百分百不是素數   

但是 此時他是素數的概率很小  很小很小  比你買的新電腦剛開啟突然爆炸了的概率還小 

所以 你可以認為他就不是個素數

那麼多嘗試幾個k  其實就是多嘗試幾個a     k值取決於 a 的值 

所以 多試幾個a   這個演算法的錯誤概率就可以降低到忽略不計

一般試個十幾次  二三十次就差不多

而且據說按照下表取a的值  在一定範圍內 可以使正確率達到百分之百   變成一個確定的演算法

  1.        if n < 1 373 653   a = 2 and 3.
  2.        if n < 9 080 191   a = 31 and 73.
  3.   if n < 4 759 123 141   a = 2, 7, and 61. 
  4.   if n < 2 152 302 898 747      a = 2, 3, 5, 7, and 11.

 

程式碼

標有註釋

模板題目51nod1106

#include<stdio.h>
#include<iostream>
#include<algorithm>
using namespace std;
typedef long long ll;
ll mul(ll a,ll b,ll mod){//高精度
    a%=mod;
    b%=mod;
    ll c=(long double)a*b/mod;
    ll ans=a*b-c*mod;
    return (ans%mod+mod)%mod;
}
ll pow_mod(ll x,ll n,ll mod){//快速冪
    ll res=1;
    while(n){
        if(n&1)
        res=mul(res,x,mod);
        x=mul(x,x,mod);
        n>>=1;
    }
    return (res+mod)%mod;
}
bool Miller_Rabbin(ll a,ll n){
    
    //把n-1  轉化成 (2^r)*d
    ll s=n-1,r=0;
    while((s&1)==0){
        s>>=1;r++;
    }
    
    //算出 2^d  存在 k 裡
    ll k=pow_mod(a,s,n);
    
    //二次探測  看變化過程中是不是等於1 或 n-1
    if(k==1)return true;
    for(int i=0;i<r;i++,k=k*k%n){
        if(k==n-1)return true;
    }
    return false;
}
bool isprime(ll n){
    //這裡可以隨機取a值進行探測  探測次數可以自己定
    //我寫了幾個我喜歡用的探測資料
    ll times=7;
    ll prime[100]={2,3,5,7,11,233,331};
    for(int i=0;i<times;i++){
        if(n==prime[i])return true;
        if(Miller_Rabbin(prime[i],n)==false)return false;//未通過探測 返回假
    }
    return true;//所有探測結束 返回真
}

 

偽素數表

常用的

 18位素數:154590409516822759
 19位素數:2305843009213693951 (梅森素數)
 19位素數:4384957924686954497

Carmichael-list:

561
1105
1729
2465
2821
6601
8911
10585
15841
29341
41041
46657
52633
62745
63973
75361
101101
115921
126217
162401
172081
188461
252601
278545
294409
314821
334153
340561
399001
410041
449065
488881
512461
530881
552721
656601
658801
670033
748657
825265
838201
852841
997633
1024651
1033669
1050985
1082809
1152271
1193221
1461241
1569457
1615681
1773289
1857241
1909001
2100901
2113921
2433601
2455921
2508013
2531845
2628073
2704801
3057601
3146221
3224065
3581761
3664585
3828001
4335241
4463641
4767841
4903921
4909177
5031181
5049001
5148001
5310721
5444489
5481451
5632705
5968873
6049681
6054985
6189121
6313681
6733693
6840001
6868261
7207201
7519441
7995169
8134561
8341201
8355841
8719309
8719921
8830801
8927101
9439201
9494101
9582145
9585541
9613297
9890881
10024561
10267951
10402561
10606681
10837321
10877581
11119105
11205601
11921001
11972017
12261061
12262321
12490201
12945745
13187665
13696033
13992265
14469841
14676481
14913991
15247621
15403285
15829633
15888313
16046641
16778881
17098369
17236801
17316001
17586361
17812081
18162001
18307381
18900973
19384289
19683001
20964961
21584305
22665505
23382529
25603201
26280073
26474581
26719701
26921089
26932081
27062101
27336673
27402481
28787185
29020321
29111881
31146661
31405501
31692805
32914441
33596641
34196401
34657141
34901461
35571601
35703361
36121345
36765901
37167361
37280881
37354465
37964809
38151361
38624041
38637361
39353665
40280065
40430401
40622401
40917241
41298985
41341321
41471521
42490801
43286881
43331401
43584481
43620409
44238481
45318561
45877861
45890209
46483633
47006785
48321001
49333201
50201089
53245921
54767881
55462177
56052361
58489201
60112885
60957361
62756641
64377991
64774081
65037817
65241793
67371265
67653433
67902031
67994641
68154001
69331969
70561921
72108421
72286501
74165065
75151441
75681541
75765313
76595761
77826001
78091201
78120001
79411201
79624621
80282161
80927821
81638401
81926461
82929001
83099521
83966401
84311569
84350561
84417985
87318001
88689601
90698401
92625121
93030145
93614521
93869665
94536001
96895441
99036001
99830641
99861985
100427041
101649241
101957401
102090781
104404861
104569501
104852881
105117481
105309289
105869401
107714881
109393201
109577161
111291181
114910489
115039081
115542505
116682721
118901521
119327041
120981601
121247281
122785741
124630273
127664461
128697361
129255841
129762001
130032865
130497361
132511681
133205761
133344793
133800661
134809921
134857801
135556345
136625941
139592101
139952671
140241361
144218341
145124785
146843929
150846961
151530401
151813201
153927961
157731841
158404141
158864833
159492061
161035057
161242705
161913961
163954561
167979421
168659569
169057801
169570801
170947105
171679561
172290241
172430401
172947529
173085121
174352641
175997185
176659201
178451857
178482151
178837201
180115489
181154701
182356993
184353001
186393481
186782401
188516329
188689501
189941761
193910977
194120389
194675041
196358977
200753281
206955841
208969201
212027401
214850881
214852609
216821881
221884001
226509361
227752993
228842209
230630401
230996949
231194965
237597361
238244041
238527745
241242001
242641153
246446929
247095361
250200721
252141121
255160621
256828321
257495641
258634741
266003101
270857521
271481329
271794601
273769921
274569601
275283401
277241401
278152381
279377281
280067761
288120421
289860481
291848401
292244833
292776121
295643089
295826581
296559361
299736181
300614161
301704985
302751505
306871201
311388337
321197185
321602401
328573477
329769721
333065305
333229141
338740417
334783585
346808881
348612265
354938221
357380101
358940737
360067201
362569201
364590721
366532321
366652201
367804801
367939585
368113411
382304161
382536001
390489121
392099401
393513121
393716701
395044651
395136505
399906001
403043257
405739681
413058601
413138881
416964241
419520241
426821473
429553345
434330401
434932961
438359041
440306461
455106601
458368201
461502097
461854261
462199681
471905281
471441001
473847121
477726145
481239361
483006889
484662529
490099681
490503601
492559141
503758801
507726901
510825601
511338241
516684961
517937581
518117041
518706721
527761081
529782121
530443201
532758241
540066241
542497201
544101481
545363281
547652161
548871961
549333121
549538081
551672221
552894301
555465601
556199281
556450777
557160241
558977761
561777121
564651361
568227241
569332177
573896881
577240273
579606301
580565233
590754385
595405201
597717121
600892993
602074585
602426161
606057985
609865201
612816751
616463809
620169409
625060801
625482001
629692801
631071001
633639097
652969351
656187001
662086041
683032801
683379841
686059921
689880801
697906561
702683101
703995733
704934361
710382401
710541481
711374401
713588401
717164449
727083001
739444021
743404663
744866305
752102401
759472561
765245881
771043201
775368901
775866001
776176261
784966297
790020001
790623289
794937601
798770161
804978721
809702401
809883361
814056001
822531841
824389441
829678141
833608321
834244501
839275921
841340521
843704401
847491361
849064321
851703301
851934601
852729121
855734401
863984881
867800701
876850801
882796321
885336481
888700681
897880321
902645857
914801665
918661501
928482241
931694401
934784929
935794081
939947009
940123801
941056273
954732853
955134181
957044881
958735681
958762729
958970545
962442001
962500561
963168193
968553181
975303121
977892241
981567505
981789337
985052881
990893569
993420289
993905641
1001152801
1027334881
1030401901
1031750401
1035608041
1038165961
1055384929
1070659201
1072570801
1093916341
1100674561
1103145121
1125038377
1131222841
1136739745
1177195201
1180398961
1189238401
1190790721
1193229577
1198650961
1200456577
1200778753
1207252621
1213619761
1216631521
1223475841
1227220801
1227280681
1232469001
1251295501
1251992281
1256855041
1257102001
1260332137
1264145401
1268604001
1269295201
1295577361
1299963601
1309440001
1312114945
1312332001
1316958721
1317828601
1318126321
1321983937
1332521065
1337805505
1348964401
1349671681
1376844481
1378483393
1382114881
1384157161
1394746081
1394942473
1404111241
1407548341
1422477001
1428966001
1439328001
1439492041
1441316269
1442761201
1490078305
1504651681
1507746241
1515785041
1520467201
1528936501
1540454761
1574601601
1576826161
1583582113
1588247851
1597821121
1626167341
1632785701
1646426881
1648076041
1659935761
1672719217
1676203201
1685266561
1688214529
1689411601
1690230241
1699279441
1701016801
1708549501
1726372441
1746692641
1750412161
1760460481
1772267281
1776450565
1778382541
1785507361
1795216501
1801558201
1803278401
1817067169
1825568641
1828377001
1831048561
1833328621
1841034961
1846817281
1848681121
1849811041
1879480513
1894344001
1899525601
1913016001
1918052065
1942608529
1943951041
1949646601
1950276565
1954174465
1955324449
1958102641
1976295241
1984089601
1988071801
2000436751
2023528501
2049293401
2064236401
2064373921
2067887557
2073560401
2080544005
2097317377
2101170097
2105594401
2107535221
2126689501
2140538401
2140699681
2301745249
2323147201
2436691321
2438403661
2444950561
2456536681
2529410281
2560600351
2598933481
2690867401
3264820001
3313196881
3480174001
3504570301
3713287801
3787491457
3835537861
4034969401
4215885697
4464573049
4562359201
4765950001
5255104513
5278692481
5781222721
5959748521
6047866621
6630702121
6916775113
7045248121
7629221377
8044161481
8251854001
8346731851
8652633601
8714965001
8976678481
9030158341
9086767201
9139810681
9293756581
9624742921
9701285761
9789244921
10119879001
10277275681
11346205609
12121569601
12173703001
12456671569
13079177569
13691148241
13946829751
14313548881
15906120889
16157879263
16765869121
17930792881
18151032901
18500666251
19742849041
21515221081
21796387201
22062297601
23224518901
23707055737
24285378001
25509692041
25749237001
26624905201
27278026129
30923424001
31876135201
34153717249
45983665729
48874811311
51436355851
51476019409
52756389001
57274147841
58094662081
59512667761
63593140801
64188507241
65700513721
67495942201
71171308081
82380774001
83946864769
100264053529
110296864801
168003672409
172018713961
173032371289
192739365541
225593397919
461574735553
464052305161
2199733160881
10028704049893
84154807001953
197531244744661
973694665856161
1436697831295441
1493812621027441
2094319836529921
2842648863161185
5778659093725441
6665161459969441
8015398984895401
8084842432668001
8188730132744161
12300849473799601
16166305446157945
32769125985828961
36629326277622001
39108343499765281
34795784213714161
37244219285276641
51116737346161201
61754693922936001
74538625278452401
91010482208190721
103183537035680641
126708084584398801
166857847951798081
186551892579535561
190232114399046721
194379756103868401
314610088213970641
335642734654849441
346413738355448401
556237362582392401
790689421836863641
1222628719906994401
1228155123631509217
1719756626091706801
2448237906710996401
2851896013395343201
3589102252820237401
4954039956700380001