51nod 1186 质数检测 V2

来源：互联网发布：陈志武的课怎么样知乎编辑：程序博客网时间：2024/05/29 17:55
#include <bits/stdc++.h>using namespace std;const int m10=1e9;const int z10=9;const int MAXL=4;struct bnum{int data[MAXL];void read(){memset(data,0,sizeof(data));char buf[32];scanf(" %s",buf);int i,j,cnt,len=strlen(buf);for(i=len-1,j=1,cnt=0;i>=0;i--,j*=10){if(j==m10){j=1;cnt++;}data[cnt]+=j*(buf[i]-'0');}}bool operator ==(const bnum &x){for(int i=0;i<MAXL;i++){if(data[i]!=x.data[i])return 0;}return 1;}bnum &operator =(const int x){memset(data,0,sizeof(data));data[0]=x;return *this;}bnum operator +(const bnum &x){bnum ret;int i,laz=0;for(i=0;i<MAXL;i++){ret.data[i]=data[i]+x.data[i]+laz;laz=ret.data[i]/m10;ret.data[i]%=m10;}return ret;}bnum operator -(const bnum &x){int i,laz=0;bnum ret;for(i=0;i<MAXL;i++){ret.data[i]=data[i]-x.data[i]-laz;if(ret.data[i]<0){ret.data[i]+=m10;laz=1;}elselaz=0;}return ret;}// assume *this<x*2bnum operator %(const bnum &x){int i;for(i=MAXL-1;i>=0;i--){if(data[i]<x.data[i]){return *this;}else if(data[i]>x.data[i]){break;}}return ((*this)-x);}bnum &div2(){int i,laz=0,tmp;for(i=MAXL-1;i>=0;i--){tmp=data[i]&1;data[i]=(data[i]+laz)>>1;laz=tmp*m10;}return *this;}bool is_odd(){return data[0]&1;}bool is_zero(){for(int i=0;i<MAXL;i++){if(data[i])return 0;}return 1;}};bnum multi(bnum &a0,bnum &b0,bnum &p){bnum tmp=a0,b=b0,ret;ret=0;while(!b.is_zero()){if(b.is_odd())ret=(ret+tmp)%p;tmp=(tmp+tmp)%p;b.div2();}return ret;}bnum powmod(bnum &a0,bnum &b0,bnum &p){bnum tmp=a0,b=b0,ret;ret=1;while(!b.is_zero()){if(b.is_odd())ret=multi(ret,tmp,p);tmp=multi(tmp,tmp,p);b.div2();}return ret;}bool miller_rabin_test(bnum &p,int iter)//检测p是否满足费马小定理，检测inter次 {int i,j,small=0,d=0;for(i=1;i<MAXL;i++){if(p.data[i])break;}if(i==MAXL){if(p.data[0]<2)return 0;if(p.data[0]==2)return 1;small=1;}if(!p.is_odd())return 0;bnum a,s,m,one,pdl;one=1;s=pdl=p-one;while(!s.is_odd())//p-1=2^d * s，计算d和奇数s{d++;s.div2();}for(i=0;i<iter;i++){//建立随机数 a，2≤a ≤p-2a=rand();if(small){a.data[0]=a.data[0]%(p.data[0]-1)+1;}else{a.data[1]=a.data[0]/m10;a.data[0]=a.data[0]%m10;}if(a==one)continue;//m=a^s %pm=powmod(a,s,p);//若 a^(s*(2^j))%p==1，则p为合数for(j=0;j<d&&!(m==one)&&!(m==pdl);j++){m=multi(m,m,p);}if(!(m==pdl)&&j>0)return 0;}return 1;}int main(){bnum x;x.read();if(miller_rabin_test(x,5))printf("Yes\n");elseprintf("No\n");}
0 0