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),因此增大检验次数就能把错误率将到很低.当然得用大整数类...
  1 #include<iostream>
  2 #include<cstdio>
  3 #include<string>
  4 #include<vector>
  5 #include<iomanip>
  6 using namespace std;
  7 const int BASE=100000000,WIDTH=8;
  8 typedef long long ll;
  9 struct BigInt{
 10     bool neg;
 11     vector<int> s;
 12     BigInt operator =(string &str){
 13         s.clear();
 14         neg=false;
 15         if(str[0]=='-'){
 16             neg=true;
 17             str[0]='0';
 18         }
 19         int x,len=(str.length()-1)/WIDTH+1;
 20         for(int i=0;i<len;i++){
 21             int end=str.length()-i*WIDTH;
 22             int start=max(0,end-WIDTH);
 23             sscanf(str.substr(start,end-start).c_str(),"%d",&x);
 24             s.push_back(x);
 25         }
 26         return *this;
 27     }
 28     BigInt operator =(int n){
 29         s.clear();
 30         if(n<0){
 31             neg=true;
 32             s.push_back(-n);
 33         }else{
 34             neg=false;
 35             s.push_back(n);
 36         }
 37         return *this;
 38     }
 39     void print(){
 40         if(neg)cout<<'-';
 41         cout<<s.back();
 42         for(int i=s.size()-2;i>=0;i--){
 43             cout<<setw(8)<<setfill('0')<<s[i];
 44         }
 45         cout<<endl;
 46     }
 47 };
 48 bool operator ==(const BigInt &a,const BigInt &b){
 49     if(a.neg!=b.neg)return false;
 50     if(a.s.size()==b.s.size()){
 51         for(int i=0;i<a.s.size();i++)
 52             if(a.s[i]!=b.s[i])return false;
 53         return true;
 54     }
 55     return false;
 56 }
 57 bool operator<(BigInt &a,BigInt &b){
 58     if(a.neg&&!b.neg)return true;
 59     if(!a.neg&&b.neg)return false;
 60     if(a.neg&&b.neg){
 61         a.neg=b.neg=false;
 62         bool ans=a<b||a==b;
 63         a.neg=b.neg=true;
 64         return !ans;
 65     }
 66     if(a.s.size()<b.s.size())return true;
 67     else if(a.s.size()>b.s.size())return false;
 68     for(int i=a.s.size()-1;i>=0;i--){
 69         if(a.s[i]==b.s[i])continue;
 70         if(a.s[i]<b.s[i])return true;
 71         return false;
 72     }
 73     return false;
 74 }
 75 bool operator >(BigInt &a,BigInt &b){return !(a==b&&a<b);}
 76 bool operator <=(BigInt &a,BigInt &b){return (a==b||a<b);}
 77 bool operator >=(BigInt &a,BigInt &b){return !(a<b);}
 78 bool operator !=(BigInt &a,BigInt &b){return !(a==b);}
 79 BigInt operator -(BigInt &a,BigInt &b);
 80 BigInt operator +(BigInt &a,BigInt &b){
 81     BigInt ans;
 82     if(a.neg&&!b.neg){
 83         a.neg=false;
 84         ans=b-a;
 85         a.neg=true;
 86         return ans;
 87     }
 88     if(!a.neg&&b.neg){
 89         b.neg=false;
 90         ans=a-b;
 91         b.neg=true;
 92         return ans;
 93     }
 94     ans.neg=a.neg;
 95     int k=0,i=0;
 96     while(i<a.s.size()||i<b.s.size()){
 97         if(i<a.s.size())k+=a.s[i];
 98         if(i<b.s.size())k+=b.s[i];
 99         ans.s.push_back(k%BASE);
100         k/=BASE;
101         i++;
102     }
103     if(k)ans.s.push_back(k);
104     return ans;
105 }
106 BigInt operator -(BigInt &a,BigInt &b){
107     BigInt ans;
108     if(!a.neg&&b.neg){
109         b.neg=false;
110         ans=a+b;
111         b.neg=true;
112         return ans;
113     }
114     if(a.neg&&!b.neg){
115         a.neg=false;
116         ans=a+b;
117         a.neg=true;
118         ans.neg=true;
119         return ans;
120     }
121     if(a.neg&&b.neg&&a>b){
122         b.neg=a.neg=false;
123         ans=b-a;
124         b.neg=a.neg=true;
125         return ans;
126     }
127     if(!a.neg&&!b.neg&&a<b){
128         ans=b-a;
129         ans.neg=true;
130         return ans;
131     }
132     ans.neg=false;
133     int k=0,i=0;
134     while(i<a.s.size()||i<b.s.size()){
135         if(i<a.s.size())k+=a.s[i];
136         if(i<b.s.size())k-=b.s[i];
137         if(k<0){
138             ans.s.push_back(k+BASE);
139             k=-1;
140         }else{
141             ans.s.push_back(k);
142             k=0;
143         }
144         i++;
145     }
146     while(ans.s.size()>1&&i>=0&&ans.s[--i]==0)ans.s.pop_back();
147     return ans;
148 }
149 BigInt operator *(BigInt &a,BigInt &b){
150     if(a.s.size()<b.s.size())return b*a;
151     BigInt ans;
152     ll k=0;
153     for(int l=0;l<a.s.size()+b.s.size()-1;l++){
154         for(int i=min(l,int(a.s.size()-1));i>=0&&l-i<b.s.size();i--){
155             k+=(ll)a.s[i]*b.s[l-i];
156         }
157         ans.s.push_back(k%BASE);
158         k/=BASE;
159     }
160     if(k)ans.s.push_back(k);
161     ans.neg=a.neg^b.neg;
162     return ans;
163 }
164 BigInt operator *(BigInt &a,int b){
165     BigInt ans;
166     ll k=0;
167     for(int l=0;l<a.s.size();l++){
168         k+=a.s[l]*b;
169         ans.s.push_back(k%BASE);
170         k/=BASE;
171     }
172     if(k)ans.s.push_back(k);
173     if((b<0&&a.neg)||(b>0&&!a.neg))ans.neg=false;
174     else ans.neg=true;
175     return ans;
176 }
177 inline void devide_2(BigInt &a){
178     int k=0;
179    // a.print();
180     for(int i=a.s.size()-1;i>=0;i--){
181         int k0=k;
182         k=a.s[i]%2==0?0:BASE/2;
183         //cout<<a.s[i]<<endl;
184         a.s[i]=k0+a.s[i]/2;
185         //cout<<a.s[i]<<endl;
186     }
187     if(a.s.back()==0)a.s.pop_back();
188 }
189 inline BigInt devide(BigInt &a,BigInt &b,BigInt &r){
190     BigInt L,R,m,t;
191     L=1;R=a;
192     if(a<b){
193         r=a;
194         return m=0;
195     }
196     while(R>L){
197         //L.print();
198         //R.print();
199         t=(R+L);
200         devide_2(t);
201         if(t==L){
202             m=t*b;
203             r=a-m;return L;
204         }
205         m=t*b;
206         //m.print();
207         if(m==a){
208             r=0;return t;
209         }else if(m<a){
210             L=t;
211         }else{
212             R=t;
213         }
214     }
215     m=t*b;
216     r=a-m;
217     return L;
218 }
219 inline BigInt operator%(BigInt &a,BigInt &b){
220     BigInt r;
221     devide(a,b,r);
222     return r;
223 }
224 inline bool Is_even(const BigInt &a){
225     if(a.s[0]%2==0)return true;
226     return false;
227 }
228 //计算n的k次方模p
229 BigInt qpow(BigInt &n,BigInt &k,BigInt &p){
230     BigInt ans,y,_n=k,a;
231     a=1;
232     //n.print();k.print();p.print();
233     ans=1;y=n;
234     while(_n>=a){
235         if(!Is_even(_n)){
236             ans=(ans*y);
237             //ans.print();
238             ans=ans%p;
239             //ans.print();
240         }
241         y=y*y;
242         //y.print();
243         y=y%p;
244         //y.print();
245         //cout<<"---------
";
246         devide_2(_n);
247        // _n.print();
248     }
249     return ans;
250 }
251 BigInt x[50];
252 int len=0;
253 void Init(){
254     int num[20];
255     for(int i=2;i<20;i++)num[i]=i;
256     for(int i=2;i<20;i++){
257         if(num[i]==0)continue;
258         x[len++]=i;
259         for(int j=2*i;j<20;j+=i)
260             num[j]=0;
261     }
262 }
263 bool Is_pri(BigInt &n){
264     if(Is_even(n)){
265         if(n.s.size()==1&&n.s[0]==2)return true;
266         return false;
267     }
268     BigInt a,b,r,u;
269     a=1;b=n-a;
270     for(int i=0;i<len;i++){
271         if(x[i]>=n)break;
272         r=b;
273         //r.print();
274         BigInt &t=x[i];
275         //t.print();
276         u=qpow(t,r,n);
277        // u.print();
278         if(u!=a)return false;
279         while(Is_even(r)){
280             devide_2(r);
281             u=qpow(t,r,n);
282             //u.print();
283             if(u!=a){
284                 if(u!=b)return false;
285                 break;
286             }
287         }
288     }
289     return true;
290 }
291 int main(){
292     //freopen("e:\out.txt","w",stdout);
293 //    BigInt a,b,r;
294 //    string s;
295 //    cin>>s;
296 //    a=s;
297 //    cin>>s;
298 //    b=s;
299 //    devide(a,b,r).print();
300 //    r.print();
301     Init();
302     BigInt n;
303     string s;
304     while(cin>>s){
305         n=s;
306         if(Is_pri(n))cout<<"Yes
";
307         else cout<<"No
";
308     }
309     return 0;
310 }