bzoj3667: RabinMiller算法

3667: Rabin-Miller算法

Input

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

Output

第一行CAS(CAS<=350，代表测试数据的组数) 
以下CAS行：每行一个数字，保证是在64位长整形范围内的正数。 
对于每组测试数据：输出Prime，代表它是质数，或者输出它最大的质因子，代表它是和数 

Miller Rabin和pollard rho模板题

卡常卡好久，最后把快速（慢速）乘改成long double才过的

#include<bits/stdc++.h>
using namespace std;
#define REP(i,st,ed) for(register int i=st,i##end=ed;i<=i##end;++i)
#define DREP(i,st,ed) for(register int i=st,i##end=ed;i>=i##end;--i)
typedef unsigned long long ll;
inline int read(){
    int x;
    char c;
    int f=1;
    while((c=getchar())!='-' && (c<'0' || c>'9'));
    if(c=='-') c=getchar(),f=-1;
    x=c^'0';
    while((c=getchar())>='0' && c<='9') x=(x<<1)+(x<<3)+(c^'0');
    return x*f;
}
inline ll readll(){
    ll x;
    char c;
    ll f=1;
    while((c=getchar())!='-' && (c<'0' || c>'9'));
    if(c=='-') c=getchar(),f=-1;
    x=c^'0';
    while((c=getchar())>='0' && c<='9') x=(x<<1ll)+(x<<3ll)+(c^'0');
    return x*f;
}
const int Times=5;
inline bool chkmax(ll &x,ll y){return (y>x)?(x=y,1):0;}
int a[4]={2,3,5,7};
ll randll(ll l,ll r){
    ll mod=(r-l+1);
    return 1ll*rand()%mod*rand()%mod+l;
}
ll work(ll x,ll y,ll mod){
    if(y>x) return x+mod-y;
    return x-y;
}
ll mul(ll x,ll y,ll mod){
    return work(x*y,(ll)((long double)x/mod*y+1e-8)*mod,mod);
}
ll ksm(ll x,ll y,ll mod){
    ll res=1;
    while(y){
        if(y&1ll) res=mul(res,x,mod);
        x=mul(x,x,mod);
        y>>=1ll;
    }
    return res;
}
bool miller_rabin(ll n){
    if(n==1 || (n>3 && n%6!=1 && n%6!=5)) return 0;
    REP(i,0,3){
        if(n==a[i]) return 1;
        if(n%a[i]==0) return 0;
    }
    ll u=n-1,cnt=0;
    while(!(u&1ll)) u>>=1ll,++cnt;
    int T=Times;
    while(T--){
        ll x=randll(2,n-1);
        ll y=ksm(x,u,n),lst;
        REP(i,1,cnt){
            lst=y,y=mul(y,y,n);
            if(y==1 && lst!=1 && lst!=n-1) return 0;
        }
        if(y!=1) return 0;
    }
    return 1;
}
ll Max;
ll gcd(ll x,ll y){
    if(!y) return x;
    return gcd(y,x%y);
}
inline ll Abs(ll x,ll y){
    if(y>x) return y-x;
    return x-y;
}
ll pollard_rho(ll n,ll c){
    ll x=rand()%n,y=x,p=1,k=2;
    for(ll i=1;p==1;++i){
        x=(mul(x,x,n)+c)%n;
        p=gcd(Abs(x,y),n);
        if(i==k) k<<=1ll,y=x;
    }
    return p;
}
void find(ll n){
    if(n==1) return;
    if(miller_rabin(n)){
        chkmax(Max,n);
        return;
    }
    ll x=n,c=0;
    while(x==n && c<=n) x=pollard_rho(n,++c);
    find(x),find(n/x);
}
int main(){
#ifndef ONLINE_JUDGE
    freopen("pollard.in","r",stdin);
    freopen("pollard.out","w",stdout);
#endif
    srand(2523452);
    int T=read();
    while(T--){
        ll x=readll();Max=0;
        find(x);
        if(Max==x) printf("Prime\n");
        else printf("%llu\n",Max);
    }
    return 0;
}