Pollard Rho 算法（分解质因数）

Pollard Rho 算法是在普通筛法上进一步优化，从 \(O(\sqrt{N})\) 变为 \(O(N^{\frac{1}{4}})\)，轻松处理 \(10^{18}\) 的数据。

渗透学习法

它叫什么：

Pollard Rho 算法。

它的地位：

大整数质因数分解的武器，经常和 Miller-Rabin 素性测试 一起使用。

它的用途：

期望时间复杂度 \(O(n^{\frac{1}{4}})\) 找到一个巨大整数 \(n\) 的一个非平凡因子（除 1、n 以外的因子）。

它的核心原理：

由于直接从 \(2\) 到 \(N\) 进行枚举时间太慢，于是我们可以随机一个数列 \(f_1,f_2,\ldots,f_N\)，我们在这个数列中找两个数 \(x,y\)。对于有可能的因子 \(p\)，若 \(|x-y|\) 是 \(p\) 的倍数，那我们就可以通过求 \(\gcd(|x-y|,N)\) 从而求出 \(p\)，然后不断的找。

具体原理：

由于生日悖论（23人中有两个人同一天生日的概率有 50%），我们通过差值求 \(p\) 的倍数是比直接随机到 \(p\) 的概率要大得多的。由于我们接下来构造的数列的特殊性，我们只需 \(O(\sqrt{p})\) 的时间就可以解决问题。

我们有一个函数 \(f(x)=(x^2+c)\ \operatorname{mod} n\)，然后我们生成一个序列：随机取一个数 \(x_1\)，令 \(x_2=f(x_1)，x_3=f(x_2),\ldots x_i=f(x_{i-1})\)，其中 \(c∈(1,n)\) 是一个随取的常数。
这个函数有一个性质，如果 \(x ≡ y\ (\operatorname{mod}\ p)\)，则 \(f(x)≡f(y)\ (\operatorname{mod}\ p)\)。然后这个函数其实是个伪随机，它其实会陷入循环，期望陷入时间为 \(O(\sqrt{p})\)，只要我们观察到这样的重复 \(x_i≡x_j\ (\operatorname{mod}\ p)\)，就可以根据 \(\gcd(|x_i-x_j,N|)\) 求出 \(N\) 的一个非平凡因子。

如何实现？

首先是 Floyd 判环。

引入一个例子。假设两个人在赛跑，A 的速度快，B 的速度慢，经过一定时间后，A 一定会和 B 相遇，且相遇时 A 跑过的总距离减去 B 跑过的总距离一定是圈长的倍数。

设 \(a=f(0)\)，\(b=f(f(0))\)，每次更新 \(a=f(a)\)，\(b=f(f(b))\)，只要检查在更新过程中 \(a\) 和 \(b\) 是否相等，如果相等了，那么就出现了环。

ll Pollard_Rho(ll N) {
	if(N == 4) return 2;
	ll c = rand() % (N - 1) + 1;
	ll a = f(0, c, N);
	ll b = f(0, c, f(0, c, N));
	while(a != b) {
		ll d = gcd(abs(b - a), N);
		if(d > 1) return d;
		a = f(a, c, N); b = f(b, c, f(b, c, N));
	}
	return N;
}

Brent 判环：Floyd 的优化。

实际上，Floyd 判环算法可以有常数上的改进。Brent 判环是从 \(k=1\) 开始持续递增 \(k\)，在第 \(k\) 轮，A 不动，B 移动 \(2^k\) 步。如果这个过程中 B 遇到了 A，则说明已经得到了环，就退出。否则让 A 移到 B 的位置，然后继续下一轮。

倍增优化：
无论是 Floyd 判环还是 Brent 判环，迭代次数都是 \(O(p)\)。然后每次迭代都需要 \(\gcd\) 一次，时间是 \(\log\ n\) 的，拖慢算法运行速度，可以通过乘法累积来减少求 \(\gcd\) 的次数。

其实就是把每次的 \(|x_i-x_j|\) 累乘起来，然后再 \(\gcd\)。原理是如果 \(|x_i-x_j|\) 中有 \(p\) 的倍数的话，它整个的乘积也必然有 \(p\) 的倍数。
如果每 \(k\) 对计算一次 \(\gcd\)，则算法复杂度降低到 \(O(\sqrt{p}+k^{-1}\sqrt{p}\log N)\)。注意到 \(k\) 和 \(\log N\) 大致同阶时，时间复杂度大约是 \(O(\sqrt{p})\)。 我们通常取 \(k=127\)。

实现：

#include <bits/stdc++.h>
#define ll long long
#define i128 __int128
using namespace std;
const int N = 2e5 + 10;
ll fac[110], cnt;

ll gcd(ll a, ll b) {return b == 0 ? a : gcd(b, a % b);}
ll qpow(ll a, ll b, ll mod) {
	ll ans = 1;
	while(b) {
		if(b & 1) ans = (i128)ans * a % mod;
		a = (i128)a * a % mod; b >>= 1;
	}
	return ans;
}
bool Miller_Rabin(ll p) {
	if(p < 2) return 0;
	if(p <= 3) return 1;
	ll d = p - 1, r = 0;
	while(!(d & 1)) r++, d >>= 1;
	for(ll k = 0; k <= 9; k++) {
		ll a = rand() % (p - 2) + 2;
		ll x = qpow(a, d, p);
		if(x == 1 || x == p - 1) continue;
		for(int i = 0; i < r - 1; i++) {
			x = (i128)x * x % p;
			if(x == p - 1) break;
		}
		if(x != p - 1) return 0;
	}
	return 1;
} // Miller Rabin 素数判定
ll Pollard_Rho(ll x) {
	ll a = 0, b = 0;
	ll c = (ll)rand() % (x - 1) + 1;
	ll val = 1;
	for(int bas = 1; ; bas <<= 1, a = b, val = 1) {
		for(int i = 1; i <= bas; i++) {
			b = ((i128)b * b + c) % x;
			val = (i128)val * abs(a - b) % x;
			if(val == 0) return x;
			if(i % 127 == 0) {
				ll d = gcd(val, x);
				if(d > 1) return d;
			}
		}
		ll d = gcd(val, x);
		if(d > 1) return d;
	}
	return x;
} // Brent 判环 + 倍增优化
void get_fac(ll x) {
	if(x < 2) return;
	if(Miller_Rabin(x)) {
		fac[++cnt] = x; return;
	}
	ll p = x;
	while(p == x) p = Pollard_Rho(x);
	get_fac(p); get_fac(x / p);
}
int main() {
	int T; scanf("%d", &T);
	while(T--) {
		ll x; scanf("%lld", &x);
		cnt = 0; get_fac(x);
		if(cnt == 1) puts("Prime");
		else {
			ll ans = 1;
			for(int i = 1; i <= cnt; i++) ans = max(ans, fac[i]);
			printf("%lld\n", ans);
		}
	}
	return 0;
}