[hiho第92周]Miller-Rabin素性测试的c++实现

证明：

如果n是素数，整数$a$ 与$n$ 互素，即$n$ 不整除$a$ ，则${a^{n - 1}} equiv 1(mod n)$ ，如果能找到一个与$n$ 互素的整数$a$ ，是的上式不成立，则可以断定$n$ 是合数，反之则不成立，这类合数我们称之为Carmichael数。当上式成立时，称$n$ 为以$a$ 为底的伪素数。

以上测试素数的方法称为fermat测试。

Miller-Rabin素性检验是在上面的基础上加上一个二次探测定理。

强伪素数：设$n - 1 = {2^s}t$ ，$2 mid t$ ，$b$ 与$n$ 互素。若${b^t} equiv 1(mod n)$ 或存在$r$ ， $0 le r le s$ 使得${b^{{2^r}t}} equiv - 1(mod n)$ ，则称n为以b为底的强伪素数。

当$n$ 为素数时，他一定是从任何数$b$ 为基的强伪素数，以$b$为基的强伪素数一定是以$b$为基的伪素数。

二次探测定理：如果p是奇素数，则 ${x^2} equiv 1(mod p)$ 的解为$x equiv 1$ 或 $x equiv p - 1(mod p)$

如果${a^{n - 1}} equiv 1(mod n)$成立，Miller-Rabin算法不是立即找另一个$a$进行测试，而是看$n-1$ 是不是偶数。如果$n-1$ 是偶数，另$u = frac{{n - 1}}{2}$，并检查是否满足二次探测定理即${a^u} equiv 1$或${a^u} equiv n - 1(mod n)$。若不满足，则为合数。

定理：若n是奇合数，则在区间$0 < b < n$ 中，最多有25%的数$b$ ，能使$n$ 是以$b$ 为基的强伪素数。

所以，结果的正确率为$1 - frac{1}{{{4^k}}}$

复杂度：$O(Slog n)$

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cstdlib>
#include<iostream>
#include<cmath>
typedef long long ll;
using namespace std;
const int S=20;
ll mod_mul(ll a,ll b,ll p){
    ll res=0;
    a%=p,b%=p;
    while(b){
        if(b&1)res=(res+a)%p;
        a=(a<<1)%p;
        b>>=1;
    }
    return res;
}
ll mod_pow(ll x,ll n,ll p){
    ll res=1;
    while(n){
        if(n&1)res=mod_mul(res,x,p);
        x=mod_mul(x,x,p);
        n>>=1;
    }
    return res;
}

bool check(ll a,ll n,ll x,ll t){//判断是否为合数 
    ll ret=mod_pow(a,x,n);
    ll last=ret;
    for(int i=1;i<=t;i++){
        ret=mod_mul(ret,ret,n);
        if(ret==1&&last!=1&&last!=n-1)return 1;
        last=ret;
    }
    if(ret!=1) return 1;//fermat测试 
    return 0;
}

bool Miller_Rabin(ll n){
    if(n<2)return 0;
    if(n==2)return 1;
    if((n&1)==0)return 0;
    ll x=n-1,t=0;
    while((x&1)==0)x>>=1,t++;
    for(int i=0;i<S;i++){
        ll a=rand()%(n-1)+1;
        if(check(a,n,x,t))return 0;//合数
    }
    return 1;
}

int main(){
    ll t,n;
    scanf("%lld",&t);
    while(t--){
        scanf("%lld",&n);
        if(Miller_Rabin(n))printf("Yes
");
        else printf("No
");
    }
    return 0;
}