洛谷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<<'\n';continue;}v.clear();fj(n);sort(v.begin(),v.end());(v.back()==n?cout<<"Prime":cout<<v.back())<<'\n';}
}