大素数判断和素因子分解（miller-rabin，Pollard_rho算法）

  1 #include<stdio.h>
  2 #include<string.h>
  3 #include<stdlib.h>
  4 #include<time.h>
  5 #include<iostream>
  6 #include<algorithm>
  7 using namespace std;
  8 
  9 
 10 //****************************************************************
 11 // Miller_Rabin 算法进行素数测试
 12 //速度快，而且可以判断 <2^63的数
 13 //****************************************************************
 14 const int S=20;//随机算法判定次数，S越大，判错概率越小
 15 
 16 
 17 //计算 (a*b)%c.   a,b都是long long的数，直接相乘可能溢出的
 18 //  a,b,c <2^63
 19 long long mult_mod(long long a,long long b,long long c)
 20 {
 21     a%=c;
 22     b%=c;
 23     long long ret=0;
 24     while(b)
 25     {
 26         if(b&1){ret+=a;ret%=c;}
 27         a<<=1;
 28         if(a>=c)a%=c;
 29         b>>=1;
 30     }
 31     return ret;
 32 }
 33 
 34 
 35 
 36 //计算  x^n %c
 37 long long pow_mod(long long x,long long n,long long mod)//x^n%c
 38 {
 39     if(n==1)return x%mod;
 40     x%=mod;
 41     long long tmp=x;
 42     long long ret=1;
 43     while(n)
 44     {
 45         if(n&1) ret=mult_mod(ret,tmp,mod);
 46         tmp=mult_mod(tmp,tmp,mod);
 47         n>>=1;
 48     }
 49     return ret;
 50 }
 51 
 52 
 53 
 54 
 55 
 56 //以a为基,n-1=x*2^t      a^(n-1)=1(mod n)  验证n是不是合数
 57 //一定是合数返回true,不一定返回false
 58 bool check(long long a,long long n,long long x,long long t)
 59 {
 60     long long ret=pow_mod(a,x,n);
 61     long long last=ret;
 62     for(int i=1;i<=t;i++)
 63     {
 64         ret=mult_mod(ret,ret,n);
 65         if(ret==1&&last!=1&&last!=n-1) return true;//合数
 66         last=ret;
 67     }
 68     if(ret!=1) return true;
 69     return false;
 70 }
 71 
 72 // Miller_Rabin()算法素数判定
 73 //是素数返回true.(可能是伪素数，但概率极小)
 74 //合数返回false;
 75 
 76 bool Miller_Rabin(long long n)
 77 {
 78     if(n<2)return false;
 79     if(n==2)return true;
 80     if((n&1)==0) return false;//偶数
 81     long long x=n-1;
 82     long long t=0;
 83     while((x&1)==0){x>>=1;t++;}
 84     for(int i=0;i<S;i++)
 85     {
 86         long long a=rand()%(n-1)+1;//rand()需要stdlib.h头文件
 87         if(check(a,n,x,t))
 88             return false;//合数
 89     }
 90     return true;
 91 }
 92 
 93 
 94 //************************************************
 95 //pollard_rho 算法进行质因数分解
 96 //************************************************
 97 long long factor[100];//质因数分解结果（刚返回时是无序的）
 98 int tol;//质因数的个数。数组小标从0开始
 99 
100 long long gcd(long long a,long long b)
101 {
102     if(a==0)return 1;//???????
103     if(a<0) return gcd(-a,b);
104     while(b)
105     {
106         long long t=a%b;
107         a=b;
108         b=t;
109     }
110     return a;
111 }
112 
113 long long Pollard_rho(long long x,long long c)
114 {
115     long long i=1,k=2;
116     long long x0=rand()%x;
117     long long y=x0;
118     while(1)
119     {
120         i++;
121         x0=(mult_mod(x0,x0,x)+c)%x;
122         long long d=gcd(y-x0,x);
123         if(d!=1&&d!=x) return d;
124         if(y==x0) return x;
125         if(i==k){y=x0;k+=k;}
126     }
127 }
128 //对n进行素因子分解
129 void findfac(long long n)
130 {
131     if(Miller_Rabin(n))//素数
132     {
133         factor[tol++]=n;
134         return;
135     }
136     long long p=n;
137     while(p>=n)p=Pollard_rho(p,rand()%(n-1)+1);
138     findfac(p);
139     findfac(n/p);
140 }
141 
142 int main()
143 {
144     //srand(time(NULL));//需要time.h头文件//POJ上G++不能加这句话
145     long long n;
146     while(scanf("%I64d",&n)!=EOF)
147     {
148         tol=0;
149         findfac(n);
150         for(int i=0;i<tol;i++)printf("%I64d ",factor[i]);
151         printf("\n");
152         if(Miller_Rabin(n))printf("Yes\n");
153         else printf("No\n");
154     }
155     return 0;
156 }