Miller-Rabin 素性检测
部分内容摘自 题解 P4718/论 Miller-Rabin 算法的确定性化 - It's LUNATIC time!)
根据费马小定理，若 \(p\) 为素数，那么对于 \(1 \leq a < p\)，都有 \(a^{p-1} \equiv 1 \pmod p\)。
因此，若存在 \(1 \leq a < p\) 使得 \(a^{p-1} \not\equiv 1 \pmod p\)，那么 \(p\) 一定不是素数。
于是有费马素性检测：在 \([1,p)\) 中随机选取几个 \(a\)，然后使用快速幂来检测是否有 \(a^{p-1} \equiv 1 \pmod p\)。如果对于每个选取的 \(a\)，\(p\) 都通过了检测，那么费马素性检测认为 \(p\) 是素数。
但是有的合数 \(p\) 满足对于 \(1 \leq a < p\)，都有 \(a^{p-1} \equiv 1 \pmod p\)。这些数被称为卡迈克尔数。因此，我们需要检测 \(p\) 的更多性质。
二次探测定理：若 \(p\) 为奇素数，则 \(x^2 \equiv 1 \pmod p\) 的解为 \(x \equiv \pm 1 \pmod p\)。
若 \(p\) 为奇数，显然 \(p-1\) 是偶数。首先令 \(d = p-1\)，然后判断 \(a^d \bmod p\) 是否为 \(1\)，若不是则 \(p\) 为合数。接下来，若 \(d\) 为偶数则令 \(d\) 除以 \(2\)，然后重复上述判断，直到 \(d\) 为奇数或 \(a^d \not\equiv 1 \pmod p\)。
如果最后一步中 \(a^d \not \equiv \pm 1 \pmod p\)，则 \(p\) 一定为合数，否则 Miller-Rabin 素性检测认为它是素数。
事实上只需要选择前 \(12\) 个素数为底数就可以适用于 \(2^{78}\) 内的素性检测。
struct MillerRabin{
int prime[20]={2,3,5,7,11,13,17,19,23,29,31,37};
inline bool test(int v,ll p){
ll d=p-1,cur=qpow(v,d,p);
if(cur!=1) return false;
while(!(d&1)){
d>>=1,cur=qpow(v,d,p);
if(cur==p-1) return true;
else if(cur!=1) return false;
}
return true;
}
inline bool check(ll x){
if(x<=40){
for(int i=0;i<12;++i) if(x==prime[i]) return true;
return false;
}
for(int i=0;i<12;++i) if(!test(prime[i],x)) return false;
return true;
}
}T;
Pollard-Rho 算法
Pollard-Rho 算法用于在 \(O(n^{1/4})\) 的复杂度计算合数 \(n\) 的某个非平凡因子。
考虑一个随机函数 \(f(x) = (x^2+c) \bmod n\)，那么取 \(f\) 的很多项相邻作差求绝对值，很大概率能够取到一个数和 \(n\) 的 GCD 不是 \(1\)。
我们倍增一下。每次取 \([2^k,2^{k+1})\) 这一段，设段头的 \(f(x)\) 值为 \(s\)，这一段里面其它的值为很多个 \(t\)，那么每次把 \(|s - t|\) 乘起来对 \(t\) 取模，如果乘了 \(127\) 次就取一下模，然后求一遍 GCD。
反正这是玄学，它能找就当它能找吧。
# include <bits/stdc++.h>
// # define int long long
const int N=100010,INF=0x3f3f3f3f;
typedef __int128 ll;
std::mt19937 rng(time(0));
inline ll qpow(ll d,ll p,ll MOD){
ll ans=1;
while(p){
if(p&1) ans=ans*d%MOD;
d=d*d%MOD,p>>=1;
}
return ans;
}
struct MillerRabin{
int prime[20]={2,3,5,7,11,13,17,19,23,29,31,37};
inline bool test(int v,ll p){
ll d=p-1,cur=qpow(v,d,p);
if(cur!=1) return false;
while(!(d&1)){
d>>=1,cur=qpow(v,d,p);
if(cur==p-1) return true;
else if(cur!=1) return false;
}
return true;
}
inline bool check(ll x){
if(x<=40){
for(int i=0;i<12;++i) if(x==prime[i]) return true;
return false;
}
for(int i=0;i<12;++i) if(!test(prime[i],x)) return false;
return true;
}
}T;
ll res;
inline ll read(void){
ll res,f=1;
char c;
while((c=getchar())<'0'||c>'9')
if(c=='-') f=-1;
res=c-48;
while((c=getchar())>='0'&&c<='9')
res=res*10+c-48;
return res*f;
}
inline ll f(ll x,ll c,ll n){
return (x*x+c)%n;
}
inline ll myabs(ll x){
return (x>0)?x:-x;
}
inline ll gcd(ll a,ll b){
return (!b)?a:gcd(b,a%b);
}
inline ll PollardRho(ll x){ // 有的时候会有 s = t，此时成环了，不需要继续找下去
ll s=0,t=0,val=1,c=rng()%(x-1)+1,d;
for(int ga=1;;ga<<=1,s=t,val=1){
for(int cur=1;cur<=ga;++cur){
t=f(t,c,x),val=myabs(s-t)*val%x;
if(cur==127){
d=gcd(val,x);
if(d>1) return d;
}
}
d=gcd(val,x);
if(d>1) return d;
}
assert(false);
return 114514;
}
inline void print(ll x){
if(x<0) putchar('-'),x=-x;
if(x>9) print(x/10);
putchar(x%10+'0');
return;
}
inline void solve(ll x){
if(x==1) return;
if(T.check(x)){
res=std::max(res,x);
return;
}
ll cur=x;
while(cur>=x) cur=PollardRho(x);
while(x%cur==0) x/=cur;
solve(cur),solve(x);
return;
}
signed main(void){
int T=read();
while(T--){
ll n=read();
res=0,solve(n);
if(res==n) puts("Prime");
else print(res),puts("");
}
return 0;
}