大素数测试和分解质因数

Prime Test http://poj.org/problem?id=1811

  1 #include<cstdio>
  2 #include<algorithm>
  3 using namespace std;
  4 typedef __int64 LL;
  5 LL mulmod(LL a,LL b,LL c) { //ret=(a*b)%c
  6     LL ret=0;
  7     for(; b; a=(a<<1)%c,b>>=1) {
  8         if(b&1) {
  9             ret=(ret+a)%c;
 10         }
 11     }
 12     return ret;
 13 }
 14 LL powmod(LL a,LL b,LL c) { //ret=(a^b)%mod
 15     LL ret=1%c;
 16     for(; b; a=mulmod(a,a,c),b>>=1) {
 17         if(b&1) {
 18             ret=mulmod(ret,a,c);
 19         }
 20     }
 21     return ret;
 22 }
 23 bool suspect(LL a,int s,LL d,LL n) {
 24     LL x=powmod(a,d,n);
 25     if(x==1) return true;
 26     while(s--) {
 27         if(x==n-1) return true;
 28         x=mulmod(x,x,n);
 29     }
 30     return false;
 31 }
 32 const int test[]= {2,3,5,7,11,13,17,19,23,-1}; // for n < 10^16
 33 bool isprime(LL n) { //Miller-Rabin 大素数测试
 34     if(n<=1||(n>2&&(!(n&1)))) return false;
 35     LL d=n-1;
 36     int s=0;
 37     while(!(d&1)) {
 38         s++;
 39         d>>=1;
 40     }
 41     for(int i=0; test[i]<n&&~test[i]; i++) {
 42         if(!suspect(test[i],s,d,n)) return false;
 43     }
 44     return true;
 45 }
 46 LL gcd(LL a,LL b){//最大公约数
 47     return b?gcd(b,a%b):a;
 48 }
 49 LL pollard_rho(LL n,LL c){//Pollard-Rho 找到n的一个约数
 50     LL d,x=rand()%n,y=x;
 51     for(LL i=1,k=2;;i++){
 52         x=(mulmod(x,x,n)+c)%n;
 53         d=gcd(y-x,n);
 54         if(d>1&&d<n) return d;
 55         if(x==y) return n;
 56         if(i==k){
 57             y=x;
 58             k<<=1;
 59         }
 60     }
 61     return 0;
 62 }
 63 LL fac[128];//每种质因数
 64 int num[128];//每种对应的个数
 65 int flen;//质因素的种数
 66 void findfac(LL n,int k) {
 67     if(n==1) return;
 68     if(isprime(n)) {
 69         fac[flen++]=n;
 70         return ;
 71     }
 72     LL p=n;
 73     int c=k;
 74     while(p>=n) {
 75         p=pollard_rho(p,c--);
 76     }
 77     findfac(p,k);
 78     findfac(n/p,k);
 79 }
 80 void gxfindfac(LL n){//大数分解质因数
 81     flen=0;
 82     findfac(n,120);
 83     sort(fac,fac+flen);
 84     num[0]=1;
 85     int k=1;
 86     for(int i=1;i<flen;i++){
 87         if(fac[i]==fac[i-1]){
 88             num[k-1]++;
 89         }
 90         else{
 91             num[k]=1;
 92             fac[k++]=fac[i];
 93         }
 94     }
 95     flen=k;
 96 }
 97 int main(){
 98     int t;
 99     LL n;
100     while(~scanf("%d",&t)){
101         while(t--){
102             scanf("%I64d",&n);
103             if(isprime(n)){
104                 puts("Prime");
105                 continue;
106             }
107             gxfindfac(n);
108             printf("%I64d\n",fac[0]);
109         }
110     }
111     return 0;
112 }

View Code