随机素数测试（Miller_Rabin算法）和求整数素因子（Pollard_rho算法）

  1 #include<cstdlib>
  2 #include<cstdio>
  3 const int maxn=10005;
  4 const int S=20;
  5 int tot;
  6 long long n;
  7 long long factor[maxn]; 
  8 long long muti_mod(long long a,long long b,long long c)
  9 {
 10     //(a*b) mod c a,b,c<2^63
 11     a%=c;
 12     b%=c;
 13     long long ret=0;
 14     while(b)
 15     {
 16         if(b&1)
 17         {
 18             ret+=a;
 19             if(ret>=c) ret-=c;
 20         }
 21         a<<=1;
 22         if(a>=c) a-=c;
 23         b>>=1;
 24     }
 25     return ret;
 26 }
 27 long long pow_mod(long long x,long long n,long long mod)
 28 {
 29     //x^n mod c
 30     if(n==1) return x%mod;
 31     int bit[90],k=0;
 32     while(n)
 33     {
 34         bit[k++]=n&1;
 35         n>>=1;
 36     }
 37     long long ret=1;
 38     for(k=k-1;k>=0;k--)
 39     {
 40         ret=muti_mod(ret,ret,mod);
 41         if(bit[k]==1) ret=muti_mod(ret,x,mod);
 42     }
 43     return ret;
 44 }
 45 bool check(long long a,long long n,long long x,long long t)
 46 {
 47     long long ret=pow_mod(a,x,n),last=ret;
 48     for(int i=1;i<=t;i++)
 49     {
 50         ret=muti_mod(ret,ret,n);
 51         if(ret==1&&last!=1&&last!=n-1) return 1;
 52         last=ret;
 53     }
 54     if(ret!=1) return 1;
 55     return 0;
 56 }
 57 bool Miller_Rabin(long long n)
 58 {
 59     long long x=n-1,t=0;
 60     while((x&1)==0) x>>=1,t++;
 61     bool flag=1;
 62     if(t>=1&&(x&1)==1)
 63     {
 64         for(int k=0;k<S;k++)
 65         {
 66             long long a=rand()%(n-1)+1;
 67             if(check(a,n,x,t)) {flag=1;break;}
 68             flag=0;
 69         }
 70     }
 71     if(flag==0||n==2) return 0;
 72     return 1;
 73 }
 74 long long gcd(long long a,long long b)
 75 {
 76     if(a==0) return 1;
 77     if(a<0) return gcd(-a,b);
 78     while(b)
 79     {
 80         long long t=a%b;a=b;b=t;
 81     }
 82     return a;
 83 }
 84 long long Pollard_rho(long long x,long long c)
 85 {
 86     long long i=1,x0=rand()%x,y=x0,k=2;
 87     while(1)
 88     {
 89         i++;
 90         x0=(muti_mod(x0,x0,x)+c)%x;
 91         long long d=gcd(y-x0,x);
 92         if(d!=1&&d!=x) return d;
 93         if(y==x0) return x;
 94         if(i==k)
 95         {
 96             y=x0;
 97             k+=k;
 98         }
 99     }
100 }
101 void findfac(long long n)  //递归分解质因数 
102 {
103     if(!Miller_Rabin(n))
104     {
105         factor[tot++]=n;
106         return;
107     }
108     long long p=n;
109     while(p>=n) p=Pollard_rho(p,rand()%(n-1)+1);
110     findfac(p);
111     findfac(n/p);
112 }
113 int main()
114 {
115     int T;
116     scanf("%d",&T);
117     while(T--)
118     {
119         scanf("%I64d",&n);
120         if(!Miller_Rabin(n))
121         {
122             printf("Prime\n");
123             continue;
124         }
125         tot=0;
126         findfac(n);
127         long long ans=factor[0];
128         for(int i=1;i<tot;i++)
129             if(factor[i]<ans) ans=factor[i];
130         printf("%I64d\n",ans);
131     }
132     return 0;
133 }