Miller_Rabin 算法 素数测试
这是一个有趣的故事。
程序员经常碰到的一个问题就是判断一个数是不是素数（经常有些无聊的数学家给出个几百位的数然后叫大家来判断是不是素数，而且被他们找出来的数还被编了号，所以经常有些什么RSA 158被分解），因为素数在计算机构建密码安全体系中占有重要的地位。平时呢就可以记一下哪些数是些素数，但本人记忆不好，只记得10以内的素数，但又不得不记一点，以后码代码还要测试样例的，所以就记一些简单的素数，比如4567是个素数，我身份证后4位是个素数，一个MM的电话号码+6是素数（所以她是我唯一一个记得电话号码的同学），23位1是个素数...
在不知道Miller_Rabin之前，一直用以下代码判断素数
1 #include<iostream> 2 #include<cmath> 3 using namespace std; 4 typedef long long ll; 5 bool fun(ll n) 6 { 7 if(n<=1) 8 return false; 9 if(n==2) 10 return true; 11 for(int i=2;i<=sqrt(n);i++) 12 if(n%i==0) 13 return false; 14 return true; 15 } 16 int main() 17 { 18 ll n; 19 while(cin>>n) 20 { 21 if(fun(n)) 22 cout<<"YES!"<<endl; 23 else 24 cout<<"NO!"<<endl; 25 } 26 return 0; 27 }
这是最基础的判断素数了，当初的优化是将2到n换到2到sqrt（n）了，我想素数的优化判断也只能到这吧。直到有一天有题目说要判断一个18位数有多少个因子，题目简单易懂，解答却要命（很多难题的题目都是很简单的，简单到你对题目只是有感觉，却毫无思路）。当初做了个简单的计算，用最基础的素数判断方法，假如碰到一个18位的素数，我的计算机就要做10^9次运算，而按照电脑每秒百万指令的运算速度，电脑每秒如果能进行10^8次运算，那么要10秒后才知道结果，但题目要求每个案例的运算时间不得超过1000ms，就是1秒，所以我能不能过题全看出题人的仁慈了，果然，最后time error在test 13；在不断优化基础算法后，比如判断奇偶数，每次i+=2，甚至i+=3；最后提交代码，但就是一直time error一直time error，time error到怀疑人生，因为我的优化时间复杂度依然是O（logn）。
我实在想不出什么方法能不经过遍历（2-sqrt（n））就能判断素数的了；直到百度后发现一种概率算法判断素数，顿时人都不好了，果然算法这东西越学越玄。
用概率怎么判断一个数是素数呢？怎么能用概率呢？既然是概率，万一错了怎么办，交题是不是还要看自己运气，交题全靠佛系？（其实概率作用很大，在一般的密码安全体系中，误判概率小于 (1/2)就可以满足安全需求。）带着恐惧去把那几位大佬的博客去看了一下下，看完之后内心平静，毫无波动，目光呆滞。
首先是学知识（发现以前数学家真可怕，动不动就去研究一些数字），还好自己离散数学还行，他们总结了很多，自己就做个小总结吧。
知识点一：一个素数x必定满足这个公式（对于一个整数a）a^(x-1)%x=1(费马小定理);证明如下
a%x=x1; 2*a%x=x2; ... (x-1)*a%x=xn;对于右边，x1，x2，x3...xn正好是1到x-1，如果相乘即（x-1）!，所以左左相乘等于右右相乘，得(x-1)!*a^(x-1)%x=(x-1)!；
两边同时约分就会得到a^(x-1)%x=1。这个公式对于所有的素数都满足，但一些合数也满足（比如Carmichael数，但Carmichael数极少，10亿里好像只有600来个），
那些合数被称为伪素数，所以一个数，如果是素数，必定满足这公式，否则就是合数。
知识点二：一个素数x，（对于一个小于它且互质的数a）a^2%x=1%x；证明如下
因为x是素数，所以只能被分解为1*x；而a^2-1=0被分解为（a+1）（a-1）%x=0,所以a只能为1或x-1。
Miller和Rabin把这些总结了下，并算了下a判错的概率为1/4,所以呢，多进行几十次判断不就可以了？嘿嘿，果然很厉害。
算法总的不多（结合快速幂，积），但看这些公式花了很久，嘿嘿，毕竟我喜欢证明这些公式。代码如下
1 #include<iostream> 2 #include<cstdlib> 3 using namespace std; 4 typedef long long ll; 5 6 ll Mul_mod(ll a,ll b,ll c)// 快速积取模 7 { 8 ll s=0; 9 while(b) 10 { 11 if(b&1) 12 s=(s+a)%c; 13 a=(a+a)%c; 14 b>>=1; 15 } 16 return s; 17 } 18 19 ll Pow_mod(ll a,ll b,ll c)// 快速幂取模 20 { 21 ll s=1; 22 while(b) 23 { 24 if(b&1) 25 s=Mul_mod(s,a,c); 26 a=Mul_mod(a,a,c); 27 b>>=1; 28 } 29 return s; 30 } 31 32 bool Miller_Rabin(ll n) 33 { 34 ll m=n-1,t=0; 35 if(n==2) //特判2 36 return true; 37 if(n<2||!(n&1)) //提前判奇偶 38 return false; 39 while(!(m&1)) 40 { 41 t++; 42 m>>=1; 43 } //m*2^t=n-1 44 for(int i=1;i<=20;i++)//测试20次 ,出错概率为4^-20 45 { 46 ll a=rand()%(n-1)+1,b=Pow_mod(a,m,n),c;//a为测试随机数 47 for(int j=1;j<=t;j++) 48 { 49 c=Mul_mod(b,b,n); 50 if(c==1&&b!=1&&b!=n-1) 51 return false; 52 b=c; 53 } 54 if(c!=1) 55 return false; 56 } 57 return true;//测试全部通过 58 } 59 60 int main() 61 { 62 ll n; 63 while(cin>>n) 64 { 65 if(Miller_Rabin(n)) 66 cout<<"YES！"<<endl; 67 else 68 cout<<"NO!"<<endl; 69 } 70 return 0; 71 }
总的来说就是我拿一些数来测你这个数，一旦你不满足这些测试，就一定不是素数，直接退出测试，判定你是合数，但是你把我测试全过了，那我就说你是素数，毕竟你过了所有测试，我只有4*-20概率说错，怕什么，哈哈，没错，你就是素数了。
以下两个参考资料：
http://www.matrix67.com/blog/archives/234
https://blog.csdn.net/ECNU_LZJ/article/details/72675595
