Miller_Rabin--大素数判定算法

一篇很好的Miller_Rabin讲解博客 费马小定理：对于一个素数$p$，对于任意整数$a$，若$gcd(a,b)=1$，则$a^{p-1}\equiv 1(modp)$ 如果要判断p是不是素数，就每次在$[2,p-1]$中随机x，如果$x^{p-1}\not\equiv 1$则p不是素数。 但这还不够精确。 二次探测定理，对于一个奇素数$p$，如果$a^2\equiv 1(modp)$，则$a\equiv 1(modp)$或$a\equiv p-1(modp)$ 我们可以把$p-1$分解成$2^t\ast u$，然后令$x[0]=a^{u}modp$，$x[i]=x[i-1{]}^{2}modp$，则$x[t]=a^{p-1}$如果$x[i]=1$，则$x[i-1]$一定等于$1$或$p-1$，否则p不为素数。如果$x[t]\ne 1$，$p$同样不为素数。

当然2要特判。 每次判断错误率为$\frac{1}{4}$，判断$s$次后错误率为$(\frac{1}{4}{)}^s$。 时间复杂度$O(slo{g}_{3}n)$ 如果p大，则乘法要用"类快速幂"计算。

long long x;
long long f[128];
long long bigrand(long long l,long long r){
    return (((long long)rand()<<31)+rand())%(r-l+1)+l;
}
long long mul(long long a,long long b,long long p){
    long long c=0;
    while(b){
        if(b&1)c=(c+a)%p;
        a=(a+a)%p;
        b>>=1;
    }
    return c;
}
long long fpm(long long a,long long b,long long p){
    long long c=1;
    while(b){
        if(b&1)c=mul(c,a,p);
        a=mul(a,a,p);
        b>>=1;
    }
    return c;
}
long long Miller_Rabin(long long x){
    long long y=x-1,t=0;
    if(x==2)return 1;
    while(!(y&1)){y>>=1;++t;}
    up(i,1,128){
        long long k=bigrand(2,x-1);
        f[0]=fpm(k,y,x);

        up(j,1,t){
            f[j]=mul(f[j-1],f[j-1],x);
            if(f[j]==1&&f[j-1]!=1&&f[j-1]!=x-1){
                return 0;
            }
        }
        if(f[t]!=1)return 0;
    }
    return 1;
}
int main(){
    long long n;
    srand(2333);
    n=read();
    while(n){
        --n;
        x=read();
        if(Miller_Rabin(x)){
            printf("Yes\n");
        }
        else printf("No\n");

    }
    return 0;
}