POJ 1811 大整数素数判断 Miller_Rabin

  1 #include <cstdio>
  2 #include <cstring>
  3 #include <cmath>
  4 #include <ctime>
  5 #include <cstdlib>
  6 #include <iostream>
  7 using namespace std;
  8 
  9 #define ll long long
 10 const int S = 10;
 11 ll ans;
 12 
 13 ll gcd(ll a,ll b){
 14     if(a < 0) return gcd(-a , b);
 15     if (b==0) return a;
 16     return gcd(b , a%b);
 17 }
 18 
 19 ll multi_mod(ll a , ll b , ll k) //a*b%k
 20 {
 21     a %= k;
 22     b %= k;
 23     ll ret = 0;
 24     //这里因为是长整数，for(i=0->b)循环一个个加太耗时，需要采用这种logn复杂度的方法
 25     while(b){
 26         if(b & 1){
 27             ret += a;
 28             ret %= k;
 29         }
 30         a <<= 1;
 31         a %= k;
 32         b >>= 1;
 33     }
 34     return ret;
 35 }
 36 
 37 ll pow_mod(ll a, ll b, ll mod) //求(a^b)%mod
 38 {
 39     if(b == 1) return a%mod;
 40     int bit[1000] , t = 0;
 41     while(b){
 42         bit[t++] = b&1;
 43         b>>=1;
 44     }
 45     //类似矩阵快速幂
 46     ll ret = 1;
 47     for(int i = t-1 ; i>=0 ; i--){
 48         ret = multi_mod(ret , ret , mod);
 49         if(bit[i]) ret = multi_mod(ret , a , mod);
 50     }
 51     return ret;
 52 }
 53 //n = x*2^t , 这里以a为底 , 检验n是否为合数
 54 bool check(ll a , ll n , ll x , ll t)
 55 {
 56     ll ret = pow_mod(a , x , n);
 57     ll last = ret; //last记录上一次的数据，保证最后结果模1时，检查到last不为n-1或1，就表示为合数
 58     for(int i = 1 ; i<=t ; i++){
 59         ret = multi_mod(ret,ret,n);
 60         if(ret == 1 && last!=1 && last!=n-1) return true;
 61         last = ret;
 62     }
 63     //始终无法找到余数为1的结果，则表示n为合数
 64     if(ret!=1) return true;
 65     return false;
 66 }
 67 //n为合数返回true
 68 bool miller_rabin(ll n)
 69 {
 70     ll x = n-1 , t = 0;
 71     while((x&1) == 0){
 72         x>>=1;
 73         t++;
 74     }
 75     bool flag = true;
 76 
 77     //保证最后得到的 x^2 mod p = 1 只有x=1或x=p-1两个解，才满足二次探测，否则有其他解就直接证明不是素数
 78     if(t >= 1 && (x&1)){
 79         //miller_robin是一种取随机测试数据的算法，这里S=20,给定20组测试数据，当然数据组数越多，正确率越大
 80         for(int i = 0 ; i<S ; i++){
 81             ll a = rand()%(n-1)+1;
 82             //二次探测在check中检测
 83             if(check(a,n,x,t)){
 84                 flag = true;
 85                 break;
 86             }
 87             flag = false;
 88         }
 89     }
 90     if(!flag || n==2)
 91         return false;
 92     else return true;
 93 }
 94 
 95 //长整数利用随机数找到其中一个因子
 96 ll Pollard_rho(ll x,ll c){
 97     ll i=1,x0=rand()%x,y=x0,k=2;
 98     while (1){
 99         i++;
100         x0=(multi_mod(x0,x0,x)+c)%x;
101         ll d=gcd(y-x0,x);
102         if (d!=1&& d!=x){
103             return d;
104         }
105         if (y==x0) return x;
106         if (i==k){
107             y=x0;
108             k+=k;
109         }
110     }
111 }
112 
113 //不断利用Pollard_rho递归来查找到n的所有素数因子
114 void find_factor(ll n)
115 {
116     if(!miller_rabin(n)){
117         ans = min(ans , n);
118         return;
119     }
120     ll p = n;
121     while(p >= n)
122         p = Pollard_rho(p , rand() % (n-1) + 1);
123     find_factor(p);
124     find_factor(n / p);
125 }
126 
127 int main()
128 {
129     srand(time(NULL));
130     int T;
131     scanf("%d" , &T);
132     while(T--){
133         ll n;
134         scanf("%lld" , &n);
135         bool flag = miller_rabin(n);
136         if(!flag) printf("Prime\n");
137         else{
138             ans = n;
139             find_factor(n);
140 
141             printf("%lld\n" , ans);
142         }
143     }
144     return 0;
145 }

 

转载于:https://www.cnblogs.com/CSU3901130321/p/4187210.html