51nod 1186 质数检测(Miller-Rabin算法)

分析:用Miller-Rabin算法，预处理一些质数，取一个质数x，若n为质数，必然有x^(n-1)=1modn，但是反之不一定，若通过了前面的检验，将n-1分解为2^d*r,r为奇数，接下来检验x^(d-1)*r mod n,以此类推除二检验，直到mod的结果为-1或不能被二整除,若其中出现mod的结果不为±1,说明n不是质数.这是因为若n为质数,x<n,>x^2=1 mod n,当且仅当x=1或x=n-1.若n通过了所有的x检验,就断定它是质数.实际上这一结论有误差,但检验k个x后错误率至多为2^(-k),因此增大检验次数就能把错误率将到很低.当然得用大整数类...

#include<iostream>
#include<cstdio>
#include<string>
#include<vector>
#include<iomanip>
using namespace std;
const int BASE=100000000,WIDTH=8;
typedef long long ll;
struct BigInt{
    bool neg;
    vector<int> s;
    BigInt operator =(string &str){
        s.clear();
        neg=false;
        if(str[0]=='-'){
            neg=true;
            str[0]='0';
        }
        int x,len=(str.length()-1)/WIDTH+1;
        for(int i=0;i<len;i++){
            int end=str.length()-i*WIDTH;
            int start=max(0,end-WIDTH);
            sscanf(str.substr(start,end-start).c_str(),"%d",&x);
            s.push_back(x);
        }
        return *this;
    }
    BigInt operator =(int n){
        s.clear();
        if(n<0){
            neg=true;
            s.push_back(-n);
        }else{
            neg=false;
            s.push_back(n);
        }
        return *this;
    }
    void print(){
        if(neg)cout<<'-';
        cout<<s.back();
        for(int i=s.size()-2;i>=0;i--){
            cout<<setw(8)<<setfill('0')<<s[i];
        }
        cout<<endl;
    }
};
bool operator ==(const BigInt &a,const BigInt &b){
    if(a.neg!=b.neg)return false;
    if(a.s.size()==b.s.size()){
        for(int i=0;i<a.s.size();i++)
            if(a.s[i]!=b.s[i])return false;
        return true;
    }
    return false;
}
bool operator<(BigInt &a,BigInt &b){
    if(a.neg&&!b.neg)return true;
    if(!a.neg&&b.neg)return false;
    if(a.neg&&b.neg){
        a.neg=b.neg=false;
        bool ans=a<b||a==b;
        a.neg=b.neg=true;
        return !ans;
    }
    if(a.s.size()<b.s.size())return true;
    else if(a.s.size()>b.s.size())return false;
    for(int i=a.s.size()-1;i>=0;i--){
        if(a.s[i]==b.s[i])continue;
        if(a.s[i]<b.s[i])return true;
        return false;
    }
    return false;
}
bool operator >(BigInt &a,BigInt &b){return !(a==b&&a<b);}
bool operator <=(BigInt &a,BigInt &b){return (a==b||a<b);}
bool operator >=(BigInt &a,BigInt &b){return !(a<b);}
bool operator !=(BigInt &a,BigInt &b){return !(a==b);}
BigInt operator -(BigInt &a,BigInt &b);
BigInt operator +(BigInt &a,BigInt &b){
    BigInt ans;
    if(a.neg&&!b.neg){
        a.neg=false;
        ans=b-a;
        a.neg=true;
        return ans;
    }
    if(!a.neg&&b.neg){
        b.neg=false;
        ans=a-b;
        b.neg=true;
        return ans;
    }
    ans.neg=a.neg;
    int k=0,i=0;
    while(i<a.s.size()||i<b.s.size()){
        if(i<a.s.size())k+=a.s[i];
        if(i<b.s.size())k+=b.s[i];
        ans.s.push_back(k%BASE);
        k/=BASE;
        i++;
    }
    if(k)ans.s.push_back(k);
    return ans;
}
BigInt operator -(BigInt &a,BigInt &b){
    BigInt ans;
    if(!a.neg&&b.neg){
        b.neg=false;
        ans=a+b;
        b.neg=true;
        return ans;
    }
    if(a.neg&&!b.neg){
        a.neg=false;
        ans=a+b;
        a.neg=true;
        ans.neg=true;
        return ans;
    }
    if(a.neg&&b.neg&&a>b){
        b.neg=a.neg=false;
        ans=b-a;
        b.neg=a.neg=true;
        return ans;
    }
    if(!a.neg&&!b.neg&&a<b){
        ans=b-a;
        ans.neg=true;
        return ans;
    }
    ans.neg=false;
    int k=0,i=0;
    while(i<a.s.size()||i<b.s.size()){
        if(i<a.s.size())k+=a.s[i];
        if(i<b.s.size())k-=b.s[i];
        if(k<0){
            ans.s.push_back(k+BASE);
            k=-1;
        }else{
            ans.s.push_back(k);
            k=0;
        }
        i++;
    }
    while(ans.s.size()>1&&i>=0&&ans.s[--i]==0)ans.s.pop_back();
    return ans;
}
BigInt operator *(BigInt &a,BigInt &b){
    if(a.s.size()<b.s.size())return b*a;
    BigInt ans;
    ll k=0;
    for(int l=0;l<a.s.size()+b.s.size()-1;l++){
        for(int i=min(l,int(a.s.size()-1));i>=0&&l-i<b.s.size();i--){
            k+=(ll)a.s[i]*b.s[l-i];
        }
        ans.s.push_back(k%BASE);
        k/=BASE;
    }
    if(k)ans.s.push_back(k);
    ans.neg=a.neg^b.neg;
    return ans;
}
BigInt operator *(BigInt &a,int b){
    BigInt ans;
    ll k=0;
    for(int l=0;l<a.s.size();l++){
        k+=a.s[l]*b;
        ans.s.push_back(k%BASE);
        k/=BASE;
    }
    if(k)ans.s.push_back(k);
    if((b<0&&a.neg)||(b>0&&!a.neg))ans.neg=false;
    else ans.neg=true;
    return ans;
}
inline void devide_2(BigInt &a){
    int k=0;
   // a.print();
    for(int i=a.s.size()-1;i>=0;i--){
        int k0=k;
        k=a.s[i]%2==0?0:BASE/2;
        //cout<<a.s[i]<<endl;
        a.s[i]=k0+a.s[i]/2;
        //cout<<a.s[i]<<endl;
    }
    if(a.s.back()==0)a.s.pop_back();
}
inline BigInt devide(BigInt &a,BigInt &b,BigInt &r){
    BigInt L,R,m,t;
    L=1;R=a;
    if(a<b){
        r=a;
        return m=0;
    }
    while(R>L){
        //L.print();
        //R.print();
        t=(R+L);
        devide_2(t);
        if(t==L){
            m=t*b;
            r=a-m;return L;
        }
        m=t*b;
        //m.print();
        if(m==a){
            r=0;return t;
        }else if(m<a){
            L=t;
        }else{
            R=t;
        }
    }
    m=t*b;
    r=a-m;
    return L;
}
inline BigInt operator%(BigInt &a,BigInt &b){
    BigInt r;
    devide(a,b,r);
    return r;
}
inline bool Is_even(const BigInt &a){
    if(a.s[0]%2==0)return true;
    return false;
}
//计算n的k次方模p
BigInt qpow(BigInt &n,BigInt &k,BigInt &p){
    BigInt ans,y,_n=k,a;
    a=1;
    //n.print();k.print();p.print();
    ans=1;y=n;
    while(_n>=a){
        if(!Is_even(_n)){
            ans=(ans*y);
            //ans.print();
            ans=ans%p;
            //ans.print();
        }
        y=y*y;
        //y.print();
        y=y%p;
        //y.print();
        //cout<<"---------
";
        devide_2(_n);
       // _n.print();
    }
    return ans;
}
BigInt x[50];
int len=0;
void Init(){
    int num[20];
    for(int i=2;i<20;i++)num[i]=i;
    for(int i=2;i<20;i++){
        if(num[i]==0)continue;
        x[len++]=i;
        for(int j=2*i;j<20;j+=i)
            num[j]=0;
    }
}
bool Is_pri(BigInt &n){
    if(Is_even(n)){
        if(n.s.size()==1&&n.s[0]==2)return true;
        return false;
    }
    BigInt a,b,r,u;
    a=1;b=n-a;
    for(int i=0;i<len;i++){
        if(x[i]>=n)break;
        r=b;
        //r.print();
        BigInt &t=x[i];
        //t.print();
        u=qpow(t,r,n);
       // u.print();
        if(u!=a)return false;
        while(Is_even(r)){
            devide_2(r);
            u=qpow(t,r,n);
            //u.print();
            if(u!=a){
                if(u!=b)return false;
                break;
            }
        }
    }
    return true;
}
int main(){
    //freopen("e:\out.txt","w",stdout);
//    BigInt a,b,r;
//    string s;
//    cin>>s;
//    a=s;
//    cin>>s;
//    b=s;
//    devide(a,b,r).print();
//    r.print();
    Init();
    BigInt n;
    string s;
    while(cin>>s){
        n=s;
        if(Is_pri(n))cout<<"Yes
";
        else cout<<"No
";
    }
    return 0;
}