寫在前面：
記錄了個人的學習過程，同時方便復習
整理自網絡
非原創部分會標明出處
by blackgryph0n
目錄 |
Miller-Rabin算法可以在O(k log2(n))的時間內檢測一個超級大的正整數n是否是素數，k為自己設定的檢測的次數
裸的Miller-Rabin算法在驗證一個數是否為素數時出錯的可能性隨着測試次數的增加而降低(不可能降為0)
- 費馬檢測
據費馬小定理，有ap ≡ a (mod p)，p為素數
而費馬小定理的逆命題是若任意正整數n滿足an ≡ a (mod n)，則n為素數
(注意，是假命題，p是素數為xp ≡ x (mod p)的充分不必要條件)
通常地，為了方便計算，使a==2
很長一段時間里，人們都認為費馬小定理的逆命題是正確的，並且有人親自驗證了a=2, p<300的所有情況
1819年有人發現了費馬小定理逆命題的第一個反例：雖然2的340次方除以341余1，但341=11*31
后來，人們又發現了561, 645, 1105等數都表明a=2時Fermat小定理的逆命題不成立。雖然這樣的數不多，但不能忽視它們的存在
於是，人們把所有能整除2^(n-1)-1的合 數n叫做偽素數，意思就是告訴人們這個素數是假的
統計表明，在前10億個自然數中共有50847534個素數，而滿足2n-1 ≡ 1 (mod n)的合數n有5597個
如果用費馬小定理的逆命題來判斷一個正整數n是不是素數，在前10億個自然數中出錯的可能性為0.011%
這個出錯的可能性還是很高的，但是仍然可以用這個技巧來排除大量的合數，這種方法就是費馬檢測
想要提高正確率可以先查詢一下偽素數表然后再使用費馬檢測，或者尋找更加准確的方法
在計算機上，計算某次冪的模一般使用快速冪[◹]快速冪
- 二次探測定理
原先在費馬小定理逆定理上的研究只考慮了a=2的情況
但是對於式子a^(n-1) mod n，取不同的a可能導致不同的結果：一個合數可能在a=2時通過了測試，但a=3時的計算結果卻排除了素數的可能
於是，人們擴展了偽素數的定義，稱滿足 a^(n-1) mod n = 1的合數n叫做以a為底的偽素數
前10億個自然數中同時以2和3為底的偽素數只有1272個，這個數目不到僅僅以2為底的偽素數的1/4。這告訴我們如果同時驗證a=2和a=3兩種情況，出錯的概率降到了0.0025%，那么以此類推，測試的數越多，出錯的可能性越小(不能降低到0)
Miller和Rabin兩個人建立了Miller-Rabin素性測試算法，基於二次探測定理
即，對於任意素數p，以及對於模p的剩余類環{1,2,...,p-1}中的任意數x，都滿足xp ≡ x (mod p)
只要存在一個數x不滿足xn ≡ x (mod n)，那么n就絕不可能是素數
為了排除大量合數，我們需要從模p的剩余類環中選取較多的數進行測試，以增強結果的可信度
每次隨機選出對於模p的剩余系{1,2,...,p-1}中的任意某個數n
將n==p-1分解為2t∗u，u為奇數
令b[0]=au mod p
令b[i]=b[i−1]2，那么b[t]=a2^t∗u
構造出b[i]后，根據二次檢測法進行判斷，如果b[i]==1，則有b[i−1]==1或者b[i−1]==p−1
如果不滿足上述條件，則p一定為合數
如果滿足了上述條件，那么p很有可能是素數
參考https://www.cnblogs.com/Norlan/p/5350243.html
然而某些合數模p的剩余類環中，對於任意1,..,p-1中的x都滿足xp ≡ x (mod p)，這類合數稱為卡邁克爾(Carmichael)數
Carmichael數是比較少的，在1~100000000范圍內的整數中，只有255個
| In number theory, a Carmichael number is a composite number n which satisfies the modular arithmetic congruence relation bn-1 ≡ 1 (mod n) for all integers b which are relatively prime to n
They are named for Robert Carmichael. The Carmichael numbers are the subset K1 of the Knödel numbers Equivalently, a Carmichael number is a composite number bn ≡ b (mod n) for which for all integers b ——Wikipedia |
| 常見卡邁克爾數： 561 1105 1729 2465 2821 6601 8911 10585 15841 …… ——bia度百科 |
裸的Miller-Rabin算法不能夠篩除這樣的合數
上代碼：
C++：
1 #include<bits/stdc++.h> 2 3 using namespace std; 4 typedef long long ll; 5 6 int const times=10; 7 8 ll ponyFE(ll a,ll b,ll c) 9 { 10 ll ans=1; 11 a%=c; 12 while(b!=0) 13 { 14 if(b&1) ans=(ans*a)%c; 15 b>>=1; 16 a=(a*a)%c; 17 } 18 return ans; 19 } 20 21 bool millerRabin(ll x) 22 { 23 if(x==2) return 1; 24 if(!(x&1)||x==1) return 0; 25 26 bool pass; 27 ll d=x-1,m; 28 while(!(d&1)) d>>=1;ll tmp=d; 29 for(int i=1;i<=times;++i) 30 { 31 d=tmp;pass=0; 32 33 m=ponyFE(rand()%(x-2)+2,d,x); 34 if(m==1) continue; 35 else for(;d<x&&d>=0;m=(m*m)%x,d<<=1) 36 if(m==x-1){pass=1;break;} 37 38 if(!pass) return 0; 39 } 40 41 return 1; 42 } 43 44 int main(int argc,char *argv[],char *enc[]) 45 { 46 int tot=0; 47 for(ll i=1;i<=100;++i) 48 if(millerRabin(i)) printf("%d\n",i),++tot; 49 printf("tot:%d\n",tot); 50 51 return 0; 52 }
Java：
1 import java.io.*; 2 import java.util.*; 3 4 class Pony 5 { 6 static int times=10; 7 8 static long ponyFE(long a,long b,long c) 9 { 10 long ans=1; 11 a%=c; 12 while(b!=0) 13 { 14 if((b&1)==1) ans=(ans*a)%c; 15 b>>=1; 16 a=(a*a)%c; 17 } 18 return ans; 19 } 20 21 static boolean millerRabin(long x) 22 { 23 if(x==2) return true; 24 if((x&1)==0||x==1) return false; 25 26 boolean pass; 27 long d=x-1,m; 28 while((d&1)==0) d>>=1;long tmp=d; 29 for(int i=1;i<=times;++i) 30 { 31 d=tmp;pass=false; 32 33 m=ponyFE((int)(Math.random())%(x-2)+2,d,x); 34 if(m==1) continue; 35 else for(;d<x&&d>=0;m=(m*m)%x,d<<=1) 36 if(m==x-1){pass=true;break;} 37 38 if(!pass) return false; 39 } 40 41 return true; 42 } 43 44 public static void main(String[] args) throws Exception 45 { 46 Scanner cin=new Scanner(System.in); 47 48 int tot=0; 49 for(long i=1;i<=100;++i) 50 if(millerRabin(i)) 51 { 52 System.out.println(i); 53 ++tot; 54 } 55 System.out.println("tot:"+tot); 56 } 57 }