【理论】Miller Rabin 和 Pollard's rho 算法

两种算法结合后可在 \(O(n^{0.25})\) 的时间复杂度内完成对一个数的质因数分解。

Miller Rabin 算法

该算法通过随机化的方式快速判定一个数是否为质数。

考虑朴素地依靠费马小定理判断一个数 \(p\) 是否为质数，即随机几个数 \(a\in[1,p)\)，检测它们是否满足 \(a^{p-1}\equiv 1\pmod p\)。

你的想法很好，但存在一种特殊的数，可以称之为 Carmichael 数，它们虽然是合数，但是却满足对于 \(\forall a\in[1,p)\)，均满足费马小定理。

那该怎么优化呢？考虑进行二次探测。

二次探测定理：若 \(p\) 为奇质数，则 \(x^2\equiv1\pmod p\) 的解为 \(x\equiv±1\pmod p\)

根据二次探测定理，如果检验 \(a^{p-1}\) 合格，那么若 \(a^{p-1}\equiv1\pmod p\)，且 \(p-1\) 为偶数，那么我们可以继续对 \(a^{\frac{p-1}{2}}\) 进行检验（注意对它进行检验时其为 \(±1\) 均可）。

不断重复将指数除以 \(2\)，直到：

\(a^{\frac{p-1}{2^k}}\ne 1\) 或 \(\frac{p-1}{2^{k}}\) 为奇数（不能继续用二次探测定理开根）

整个过程如下：

\(1.\) 先判断 \(a^{p-1}\) 模意义下是否为 \(1\)。

\(2.\) 不断将指数除 \(2\) 直到其不能除。

\(3.\) 判断此时的指数是否为 \(±1\)。

这样对于单个 \(a\) 的判断就做完了。

我们一般来说随机 \(20\) 个（保险） \(a\)，都判断一下，如果全部通过，那么便可以将其视为质数。

注意二次探测定理只针对奇素数，所以需要特判 \(p=2\) 的情况。

这样 Miller-Rabin 就做完了。

代码

bool mr(ll x){
	if(x<=1) return 0;
	if(x<=3) return 1; 
	rep(i,1,20){
		ll c=rand()%(x-3)+2;
		ll now=x-1;
		ll qwq=qp(c,now,x);
		if(qwq!=1) return 0;
		while(now%2==0){
			if(qp(c,now,x)!=1) break;
			now/=2;
		}
		ll tmp=qp(c,now,x);
		if(tmp!=1&&tmp!=x-1) return 0; 
	}	
	return 1;
}

Pollard's rho 算法

以下简称 Pr 算法。

Pr 算法所涉及的一个重要的结论就是生日悖论。

一个班有 \(n\) 个人，存在至少两人生日重复的概率是多少？

结果有点出乎我们的意料，对于 \(n=23\)， 这个概率接近 \(50\%\)，对于 \(n=60\)，这个概率约为 \(0.9999\%\)，几乎可以和 \(100\%\) 划上等号。

换个角度，要想存在至少两人生日重复（设一年有 \(n\) 天），那么期望班里有 \(O(\sqrt n)\) 个人时，出现重复。

我们可以根据这个设计出分解质因数的方法。

随机生成一个数列（所有项 \(\bmod n\)），那么这个数列在模 \(m\) 的剩余系下期望 \(O(\sqrt m)\) 次出现重复。

数列的生成方式是随机 \(x_0,c\)，之后令 \(x_i=(x_{i-1}^2+c)\bmod n\)。

设 \(n\) 的最小质因子为 \(p\)，那么生成的数列将在模 \(p\) 的剩余系下期望 \(O(\sqrt p)\) 次出现重复。

则我们为这个生成的数列中添加一项之后，我们都讲这个新的数与前面的数分别相减，并与 \(n\) 取一次 \(\gcd\)，若结果不为 \(0 or 1\)，那么就成功地将 \(n\) 分解了。

但如果我们生成了这个数列的很多项之后仍没有找到重复怎么办？

只有一种可能，就是生成的数列进入了死循环（形状类似希腊字母 rho ρ）。

如何判死循环？

需要用到一种名为 Floyd 的判环方法，考虑每次取生成的数列 \(b_k\) 和 \(b_{2k}\)。若这两者在模 \(n\) 的剩余系下相同，那么可以得到它们之间跨过了环长度的整数倍，则出现了死循环，可以直接终止。

所以我们维护两个变量 \(X,Y\)，\(X\) 每次走一步，\(Y\) 每次走两步。这样的话，我们就实现了判环，同时，我们还可以优化找两数相减的过程，将任意两对数的相减转化为每次 \(X,Y\) 的相减，虽然这回略过一些可能正确的答案，但实测效率更高。

这样，我们可以通过 Pr 算法，通过极大概率（如果找不到就再随机找一遍）找到 \(n\) 的一个因子（注意不是质因子，其内部可能有很多质因子）\(q\)，那么接下来我们需要做的就是 \(Pr(n/q)\) 和 \(Pr(q)\)，直到完全分解为止。

这样就做完了，但这样并不足以通过模板题，所以还需要在 Pr 算法内部进行优化。

我们考虑每次找到两个项的差之后都需要做一遍 \(\gcd\)，很慢，思考一下会发现我们可以把几次作差的结果乘起来再一起做 \(\gcd\)，那这个取 \(\gcd\) 的间隔是多久呢？考虑倍增，间隔一开始为 \(1\)，之后不断扩大，最后扩大到一个阈值 \(\omega\) 时停止扩大，间隔不再扩大。

一般取 \(\omega=64\) 或 \(\omega=128\)。

这样的话就能过通过本题了。

代码

#include<bits/stdc++.h>
using namespace std;
#define rep(i,a,b) for(register int i=a,i##end=b;i<=i##end;++i)
typedef long long ll;
ll mul(ll a,ll b,ll mod){return (__int128)a*b%mod;}
ll gcd(ll x,ll y){return !y?x:gcd(y,x%y);}
ll qp(ll b,ll p,ll mod){
	ll ans=1,base=b%mod;
	while(p){
		if(p&1) ans=mul(ans,base,mod);
		base=mul(base,base,mod);
		p>>=1;
	}
	return ans;
}
bool mr(ll x){
	if(x<=1) return 0;
	if(x<=3) return 1; 
	rep(i,1,20){
		ll c=rand()%(x-3)+2;
		ll now=x-1;
		ll qwq=qp(c,now,x);
		if(qwq!=1) return 0;
		while(now%2==0){
			if(qp(c,now,x)!=1) break;
			now/=2;
		}
		ll tmp=qp(c,now,x);
		if(tmp!=1&&tmp!=x-1) return 0; 
	}	
	return 1;
}
ll to(ll x,ll mod,ll c){return (mul(x,x,mod)+c)%mod;}
ll pr(ll x){
	ll c=rand()%x+1;
	ll X=rand()%x+1,Y=to(X,x,c),cnt=0,prod=1,step=1;
	while(X!=Y){
		prod=mul(prod,abs(X-Y),x);
		++cnt;
		if(cnt==step){
			cnt=0;
			if(gcd(prod,x)!=1) return gcd(prod,x);
			prod=1;
			step<<=1;
			if(step>64) step=64;
		}
		X=to(X,x,c);
		Y=to(to(Y,x,c),x,c);
	}
	return gcd(prod,x);
}
ll mx;
void getp(ll x){
	if(x==1) return;
	if(mr(x)){
		mx=max(mx,x);
		return;
	}
	ll tmp=pr(x);
	getp(tmp);
	getp(x/tmp);
}
ll t,n;
int main(){
	cin>>t;
	while(t--){
		cin>>n;
		if(mr(n)) puts("Prime");
		else{
			mx=0;
			getp(n);
			printf("%lld\n",mx);
		}
	}
}