洛谷P4718 【模板】Pollard-Rho算法

虽然很久以前就听说过PR算法，但前几天第一次打。

首先miller rabin判断素数，不在复杂度瓶颈。
pollard rho倍增环长，复杂度是(O(n^{frac{1}{4}} log n))的。
然而这样复杂度较高，比较难过加强后的数据。

可以考虑每次倍增时把乘积存下来，然后在增加环长时一次gcd，这样应该是(O(n^{frac{1}{4}}))

#pragma GCC optimize(3)
#pragma GCC optimize("Ofast")
#pragma GCC optimize("inline")

#include <bits/stdc++.h>
using namespace std;
typedef long long i64;
typedef unsigned long long u64;
typedef vector<i64> vi64;
i64 fast_pow(i64 x,i64 y,i64 z){
    //cerr<<"fp"<<x<<" "<<y<<" "<<z<<endl;
    i64 ret=1;
    for (; y; y>>=1,x=(__int128)x*x%z)
        if (y&1) ret=(__int128)ret*x%z;
    return ret;
}
u64 RAND(){
    return ((u64)rand()*rand()+rand());
}
const int sp[]={2,3,5,7,11,13,17,19,23,29,31,37,43,47};
bool is_prime(i64 n){
    //cerr<<"is_prime"<<n<<endl;
    if (n==2) return 1;
    if (~n&1) return 0;
    i64 u=n-1;
    int b=0;
    while (~u&1){
        u>>=1;
        ++b;
    }
    for (int i=0; i<12&&sp[i]<n; ++i){
        i64 now=fast_pow(sp[i],u,n),la=now;
        //cerr<<a<<" "<<u<<" "<<now<<" "<<i<<endl;
        if (now!=1)
        for (int i=1; i<=b; ++i){
            now=(__int128)now*now%n;
            if (now==1){
                if (la!=n-1) return 0;
                break;
            }
            la=now;
        }
        //cerr<<"now"<<now<<endl;
        if (now!=1) return 0;
    }
    return 1;
}
i64 fix(i64 a,i64 n){
    return a<0?a+n:a;
}
i64 gcd(i64 a,i64 b){
    int shift=__builtin_ctzll(a|b);
    b>>=__builtin_ctzll(b);
    while (a){
        a>>=__builtin_ctzll(a);
        if (a<b) swap(a,b);
        a-=b;
    }
    return b<<shift;
}
i64 pollard_rho(i64 n,i64 c){
    //cerr<<"pollard_rho"<<n<<" "<<c<<endl;
    i64 i=1,k=2,x=RAND()%(n-1)+1,y=x,g=1;
    while (1){
        ++i;
        x=((__int128)x*x+c)%n;
        //i64 tmp=gcd(fix(x-y,n),n);
        //cerr<<x<<" "<<y<<endl;
        //if (tmp>1) return tmp;
        g=(__int128)fix(x-y,n)*g%n;
        if (i==k){
        	g=gcd(g,n);
        	if (g>1){
        		if (g==n){
        			x=y;
        			while (g==1){
        				x=((__int128)x*x+c)%n;
        				g=gcd(fix(x-y,n),n);
        			}
        		}
        		return g;
        	}
            y=x;
            k<<=1;
        }
    }
}
vi64 v;
void fj(i64 n){
    if (n==1) return;
    if (is_prime(n)){
        //cerr<<"N"<<n<<endl;
        v.push_back(n);
        return;
    }
    i64 p=n;
    while (p>=n){
        p=pollard_rho(n,RAND()%n);
    }
    fj(p);
    fj(n/p);
}
typedef vector<pair<i64,int> > vpi64i;
vpi64i vv;
int k,p;
int fp(int x,int y){
    int ret=1;
    for (; y; y>>=1,x=(i64)x*x%p) if (y&1) ret=(i64)ret*x%p;
    return ret;
}
int C(int x){
    return ((u64)x*(x+1)>>1)%p;
}
int mul(int x,int y){
    return (i64)x*y%p;
}
int add(int x,int y){
    return (x+=y)>=p?x-p:x;
}
int main(){
    int t;
    cin>>t;
    while (t--){
        i64 n;
        cin>>n;
        //FastDiv fd(p);
        if (n==1){
            cout<<1<<'
';
            continue;
        }
        v.clear();
        fj(n);
        sort(v.begin(),v.end());
        (v.back()==n?cout<<"Prime":cout<<v.back())<<'
';
    }
}