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
2
13
134
8897
1234567654321
1000000000000
Sample Output
Prime
Prime
67
41
4649
5
Prime
67
41
4649
5
HINT
数据范围：
保证cas<=350，保证所有数字均在64位长整形范围内。
Source
题解：
这是miller rabin和pollard rho的裸题，具体见:http://blog.csdn.net/thy_asdf/article/details/51347390
code:
1 #include<cstdio> 2 #include<iostream> 3 #include<cmath> 4 #include<cstring> 5 #include<algorithm> 6 #define lowbit(x) ((x)&(-(x))) 7 using namespace std; 8 typedef long long int64; 9 char ch; 10 bool ok; 11 void read(int &x){ 12 ok=0; 13 for (ch=getchar();!isdigit(ch);ch=getchar()) if (ch=='-') ok=1; 14 for (x=0;isdigit(ch);x=x*10+ch-'0',ch=getchar()); 15 if (ok) x=-x; 16 } 17 void read(int64 &x){ 18 ok=0; 19 for (ch=getchar();!isdigit(ch);ch=getchar()) if (ch=='-') ok=1; 20 for (x=0;isdigit(ch);x=x*10+ch-'0',ch=getchar()); 21 if (ok) x=-x; 22 } 23 int cases; 24 int64 n; 25 const int prime[10]={2,3,5,7,11,13,17,19,23,29}; 26 int64 mul(int64 a,int64 b,int64 mod){ 27 int64 d=((long double)a/mod*b+1E-8); 28 int64 res=a*b-d*mod; 29 if (res<0) res+=mod; 30 return res; 31 /*int64 t; 32 for (t=0;b;b>>=1,a<<=1,a%=mod) if (b&1) t=(t+a)%mod; 33 return t;*/ 34 } 35 int64 ksm(int64 a,int64 b,int64 mod){ 36 int64 t; 37 for (t=1;b;b>>=1,a=mul(a,a,mod)) if (b&1) t=mul(t,a,mod); 38 return t; 39 } 40 bool miller_rabin(int64 n){ 41 for (int i=0;i<10;i++) if (n==prime[i]) return true; 42 if (!(n&1)) return false; 43 int64 bit=lowbit(n-1),s=0,d; 44 while (bit!=1) s++,bit>>=1; 45 d=(n-1)>>s; 46 for (int i=0;i<10;i++){ 47 int64 x=ksm(prime[i],d,n); 48 for (int j=1;j<=s;j++){ 49 int64 xx=mul(x,x,n); 50 if (xx==1&&x!=n-1&&x!=1) return false; 51 x=xx; 52 } 53 if (x!=1) return false; 54 } 55 return true; 56 } 57 int cnt; 58 int64 list[50]; 59 int64 random(int64 lim){return ((1LL*rand()<<31)+rand())%lim;} 60 int64 f(int64 x,int64 mod,int64 c){return (mul(x,x,mod)+c+mod)%mod;} 61 int64 pollard_rho(int64 n,int64 c){ 62 int64 x,y,d=1; x=random(n),y=f(x,n,c); 63 while (d==1){ 64 d=__gcd(abs(x-y),n); 65 x=f(x,n,c),y=f(f(y,n,c),n,c); 66 } 67 return d; 68 } 69 void work(int64 n){ 70 if (miller_rabin(n)){list[++cnt]=n;return;} 71 int64 d=pollard_rho(n,random(n-1)); 72 while (d==n||d==1) d=pollard_rho(n,random(n)); 73 work(d),work(n/d); 74 } 75 void decompose(int64 n){ 76 cnt=0,work(n),sort(list+1,list+cnt+1); 77 printf("%lld\n",list[cnt]); 78 } 79 int main(){ 80 srand(19990617); 81 for (read(cases);cases;cases--){ 82 read(n); 83 if (miller_rabin(n)) puts("Prime"); 84 else decompose(n); 85 } 86 return 0; 87 }
