poj 1181(大数判素数 ，分解) 

  1 //miller_rabin 判断一个大数是否是素数
  2 //pollard_rho 大数因子分解
  3 #include<cstdio>
  4 #include<cstdlib>
  5 #include<cstring>
  6 #include<cmath>
  7 #include<algorithm>
  8 
  9 using namespace std;
 10 
 11 #define LL long long
 12 const LL Max = (LL) 1 << 62;
 13 LL p[10]={2,3,5,7,11,13,17,19,23,29};
 14 
 15 //计算a*b%n
 16 inline LL multi_mod(LL a,LL b,LL mod)
 17 {
 18     LL sum=0;
 19     while(b)
 20     {
 21         if(b&1)    sum = (sum + a) % mod;
 22         a <<= 1;
 23         b >>= 1;
 24         if (a >= mod) a %= mod;
 25     }
 26     return sum;
 27 }
 28 
 29 //计算a^b%n;
 30 inline LL quick_mod(LL a,LL b,LL mod)
 31 {
 32     LL sum = 1;
 33     while(b)
 34     {
 35         if(b & 1) sum = multi_mod(sum, a, mod);
 36         a = multi_mod(a, a, mod);
 37         b >>= 1;
 38     }
 39     return sum;
 40 }
 41 
 42 bool miller_rabin(LL n)
 43 {
 44     LL u,m,buf;
 45     int k = 0;
 46     //将n分解为m*2^k
 47     if(n==2)   return true;
 48     if(n<2 || !(n&1))    return false;
 49     m = n - 1;
 50     while(!(m&1))
 51         k++,m >>= 1;
 52     for(int i = 0; i < 9; ++i)
 53     {
 54         if(p[i] >= n) return true;
 55         u = quick_mod(p[i], m, n);
 56         if(u==1) continue;
 57         for(int j = 0; j < k; ++j)
 58         {
 59             buf = multi_mod(u, u, n);
 60             if(buf == 1 && u != 1 && u != n - 1)
 61                 return false;
 62             u = buf;
 63         }
 64         //如果p[i]^(n-1)%n!=1那么n为合数
 65         if(u - 1) return false;
 66     }
 67     return true;
 68 }
 69 
 70 LL gcd(LL a, LL b)
 71 {
 72     return b == 0 ? a: gcd(b, a%b);
 73 }
 74 //寻找n的一个因子，该因子并不一定是最小的，所以下面要二分查找最小的那个因子
 75 LL pollard(LL n )
 76 {
 77    LL x,y,c=0,d,i=1,k=2;
 78    while(c==0 || c==2 ) c = abs( rand())%(n-1) + 1;
 79    x = y = (rand( )%( n-1 )) + 1;
 80    do
 81    {
 82       i++;
 83       d = gcd( y + n - x, n );
 84       if( d >1 && d < n )
 85           return d; 
 86       if( i == k )
 87       {
 88          y = x; k <<= 1;        
 89       }    
 90       x = (multi_mod( x , x, n )+n-c )%n;
 91     }while( x != y );
 92     return n; 
 93 }
 94 
 95 LL pollard_min(LL n)
 96 {
 97     LL p,a,b=Max;
 98     if(n==1)    return Max;
 99     if(miller_rabin(n))    return n;
100     p = pollard(n);
101     a = pollard_min(p);//二分查找
102     b = pollard_min(n / p);
103     return a < b ? a : b;
104 }
105 
106 int main(void)
107 {
108     LL T;
109     LL n;
110     scanf("%lld",&T);
111     while(T--)
112     {
113         scanf("%lld",&n);
114         if(miller_rabin(n))
115         {
116             puts("Prime");
117         }
118         else printf("%lld\n",pollard_min(n));
119     }
120     return 0;
121 }

 

posted @ 2013-05-28 16:29  Missa  阅读(553)  评论(0)    收藏  举报