bzoj3667: Rabin-Miller算法

Description

Input

第一行：CAS,代表数据组数（不大于350），以下CAS行，每行一个数字，保证在64位长整形范围内，并且没有负数。你需要对于每个数字：第一，检验是否是质数，是质数就输出Prime
第二，如果不是质数，输出它最大的质因子是哪个。

Output

第一行CAS(CAS<=350，代表测试数据的组数)

以下CAS行：每行一个数字，保证是在64位长整形范围内的正数。

对于每组测试数据：输出Prime，代表它是质数，或者输出它最大的质因子，代表它是和数

Sample Input

6
2
13
134
8897

1234567654321
1000000000000

Sample Output

Prime
Prime
67
41
4649

5

HINT

数据范围：

保证cas<=350，保证所有数字均在64位长整形范围内。

Source

题解：

这是miller rabin和pollard rho的裸题，具体见:http://blog.csdn.net/thy_asdf/article/details/51347390

code:

 #include<cstdio>

 #include<iostream>

 #include<cmath>

 #include<cstring>

 #include<algorithm>

 #define lowbit(x) ((x)&(-(x)))

 using namespace std;

 typedef long long int64;

 char ch;

 bool ok;

 void read(int &x){

     ok=;

     for (ch=getchar();!isdigit(ch);ch=getchar()) if (ch=='-') ok=;

     for (x=;isdigit(ch);x=x*+ch-'',ch=getchar());

     if (ok) x=-x;

 }

 void read(int64 &x){

     ok=;

     for (ch=getchar();!isdigit(ch);ch=getchar()) if (ch=='-') ok=;

     for (x=;isdigit(ch);x=x*+ch-'',ch=getchar());

     if (ok) x=-x;

 }

 int cases;

 int64 n;

 const int prime[]={,,,,,,,,,};

 int64 mul(int64 a,int64 b,int64 mod){

     int64 d=((long double)a/mod*b+1E-);

     int64 res=a*b-d*mod;

     if (res<) res+=mod;

     return res;

     /*int64 t;

     for (t=0;b;b>>=1,a<<=1,a%=mod) if (b&1) t=(t+a)%mod;

     return t;*/

 }

 int64 ksm(int64 a,int64 b,int64 mod){

     int64 t;

     for (t=;b;b>>=,a=mul(a,a,mod)) if (b&) t=mul(t,a,mod);

     return t;

 }

 bool miller_rabin(int64 n){

     for (int i=;i<;i++) if (n==prime[i]) return true;

     if (!(n&)) return false;

     int64 bit=lowbit(n-),s=,d;

     while (bit!=) s++,bit>>=;

     d=(n-)>>s;

     for (int i=;i<;i++){

         int64 x=ksm(prime[i],d,n);

         for (int j=;j<=s;j++){

             int64 xx=mul(x,x,n);

             if (xx==&&x!=n-&&x!=) return false;

             x=xx;

         }

         if (x!=) return false;

     }

     return true;

 }

 int cnt;

 int64 list[];

 int64 random(int64 lim){return ((1LL*rand()<<)+rand())%lim;}

 int64 f(int64 x,int64 mod,int64 c){return (mul(x,x,mod)+c+mod)%mod;}

 int64 pollard_rho(int64 n,int64 c){

     int64 x,y,d=; x=random(n),y=f(x,n,c);

     while (d==){

         d=__gcd(abs(x-y),n);

         x=f(x,n,c),y=f(f(y,n,c),n,c);

     }

     return d;

 }

 void work(int64 n){

     if (miller_rabin(n)){list[++cnt]=n;return;}

     int64 d=pollard_rho(n,random(n-));

     while (d==n||d==) d=pollard_rho(n,random(n));

     work(d),work(n/d);

 }

 void decompose(int64 n){

     cnt=,work(n),sort(list+,list+cnt+);

     printf("%lld\n",list[cnt]);

 }

 int main(){

     srand();

     for (read(cases);cases;cases--){

         read(n);

         if (miller_rabin(n)) puts("Prime");

         else decompose(n);

     }

     return ;

 }

秒客网