【朝夕的ACM笔记】目录与索引
Miller Rabin素数判定
前言
Miller Rabin素性检验是一种素数判定的法则，由CMU的教授Miller首次提出，并由希大的Rabin教授作出修改，变成了今天竞赛人广泛使用的一种算法，故称Miller Rabin素性检验。
该算法本质上是一种随机化算法，能在 O(klog^3n) 的时间复杂度下快速判断出一个数是否是素数，但具有一定的错误概率。不过在一定数据范围内，通过一些技巧可以使其不出错。
Fermat素性检验
由费马小定理我们得知：对任意素数 p 和与其互质的整数 a ， a^{p-1}\equiv1(mod\ p) 。
所以我们会思考，这个定理是否可以反过来使用？
即，随机选取一小于 p 的整数 a ，若 a^{p-1}\equiv1(mod\ p) ，则是否可以证明 p 是一素数？
很遗憾，答案是不行，因为存在一些合数满足这个等式。
这些合数被称为费马伪素数，例如最小的费马伪素数为341。
那我们又考虑，如果多随机选取几个 a 的话能否提高正确率？
是的，通过多选取几个底数，我们很大程度上降低了错误的概率，比如341就成功被筛去了。
但仍旧存在极少的一些合数，即便遍历 [2,p-1] 的每一个数字作为底数，也无法筛去。
这样的合数被称为卡迈克尔数，在一亿内有255个，最小的卡迈克尔数为561。
若 n 为卡迈克尔数，则 2^n-1 也是卡迈克尔数，故其个数是无穷的。
Miller Rabin素性检验
单单费马小定理行不通，所以我们考虑引入新的定理来提高检测的正确性。
二次探测定理：对于质数 p ，若 x^2\equiv1(mod\ p) ，则小于 p 的解只有两个， x_1=1,x_2=p-1 。
定理证明：
x^2-1\equiv0(mod\ p)——(x+1)(x-1)\equiv0(mod\ p) 。
故 p|(x+1)(x-1) ，又因为 p 是质数，故要么 (x+1)(x-1)=0 ，要么 (x+1)(x-1) 是 p 的倍数，可得唯二解。
这个定理有什么用呢？
由Fermat检测得到 a^{p-1}\equiv1(mod\ p) ，且 p-1 为偶数（否则 p 为偶数直接特判筛掉），则 a^{p-1} 就相当于 x^2 。
将其拆分为 (a^{\frac{p-1}{2}})^2\equiv1 ，就可以用二次检测定理来判断了。
如果 a^{\frac{p-1}{2}} 在 mod\ p 情况下的解不是1或者 p-1 ，那 p 就不是素数。
如果 a^{\frac{p-1}{2}}\equiv1(mod\ p) ，那可以模仿之前的操作，再进行一轮检验，变成判断 a^{\frac{p-1}{4}} ，直到最后变成奇数。
也就是说，我们可以将 p-1=u\times2^t(u为奇数) ，对 a^u,a^{u\times2},a^{u\times2^2}… 这一系列数进行检验，他们的解要么全是1，要么出现 p-1 后全是1（之前不能有1），否则就不是素数。当然，要注意 p-1 不能出现在最后一个数，否则就连费马小定理都不满足了。
还要注意这过程中不能产生 p 的倍数。
总结一下，该检验的思路是：
1）先特判筛掉3以下的数与偶数。
2）将待检验数 n ， n-1 化为 u\times2^t 。
3）选取多个底数，分别对 a^u,a^{u\times2},a^{u\times2^2}… 进行检验，判断其解是否全为1，或在非最后一个数的情况下出现 p-1 。
4）如果都满足，我们就认为这个数是素数。
代码解决
那么思路的问题解决了，接下来考虑代码的实现。
首先是底数的选取：Miller Rabin检测依旧有错误的概率，虽然微乎其微，但我们显然不想接受。但通过选取适合的底数，可以避免这一情况的发生。
在int范围内，选取 \{2,7,61\} 三个数作为底数可以保证100%正确。
在long long范围内，选取 \{2,325,9375,28178,450775,9780504,1795265022\} 七个数作为底数可保证100%正确。
故正常情况下时间复杂度最多为 O(log^3n) ，常数为7。
其次，由于这个过程中可能要计算long long型变量的平方，所以要考虑数据溢出的问题，解决方案是快速乘。
这里提供一个 O(1) 的快速乘，原理是用溢出来解决溢出。
//ld=long double,ull=unsigned long long
ll qmul(ll a,ll b,ll mod)a*b%mod
{
ll c=(ld)a/mod*b;
ll res=(ull)a*b-(ull)c*mod;
return (res+mod)%mod;
}
总参考代码：
ll qmul(ll a,ll b,ll mod)//快速乘
{
ll c=(ld)a/mod*b;
ll res=(ull)a*b-(ull)c*mod;
return (res+mod)%mod;
}
ll qpow(ll a,ll n,ll mod)//快速幂
{
ll res=1;
while(n)
{
if(n&1) res=qmul(res,a,mod);
a=qmul(a,a,mod);
n>>=1;
}
return res;
}
bool MRtest(ll n)//Miller Rabin Test
{
if(n<3||n%2==0) return n==2;//特判
ll u=n-1,t=0;
while(u%2==0) u/=2,++t;
ll ud[]={2,325,9375,28178,450775,9780504,1795265022};
for(ll a:ud)
{
ll v=qpow(a,u,n);
if(v==1||v==n-1||v==0) continue;
for(int j=1;j<=t;j++)
{
v=qmul(v,v,n);
if(v==n-1&&j!=t){v=1;break;}//出现一个n-1，后面都是1，直接跳出
if(v==1) return 0;//这里代表前面没有出现n-1这个解，二次检验失败
}
if(v!=1) return 0;//Fermat检验
}
return 1;
}