Miller-Rabin算法
写在前面：
记录了个人的学习过程，同时方便复习
整理自网络
非原创部分会标明出处
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 }
