Pollard-Rho 学习笔记

距离省选只剩 2days 了，学学概念就差不多得了吧（

一言以蔽之，Pollard-Rho 是一种基于随机的、能够在期望 $\Theta(n^{0.25})$ 的时间内对一个数 $n$ 分解质因数的算法。

那么 Pollard-Rho 究竟是如何运用随机的性质实现这一功能的呢？首先考虑一前置知识，如何使用随机化的知识判定一个数 $n$ 是否是质数，也就是传说中的 Miller-Rabin 素性测试。

暴力随机 $n$ 的因数显然是不可行的，毕竟随机到的概率太小了。考虑费马小定理：对于质数 $p$，有 $a^{p-1}\equiv 1\pmod{p}$，其中 $a\not\equiv 0\pmod{p}$，因此，我们如果找到了一个 $a\not\equiv 0\pmod{p}$，满足 $a^{p-1}\not\equiv 1\pmod{p}$，那么即可说明 $p$ 不是质数。对于某些 $p$，运用此种方法确实能以较高的概率判定一个数是否是质数，因为对于固定的合数 $p$，如果至少存在一个 $a$ 满足 $a^{p-1}\not\equiv 1\pmod{p}$，那么满足 $a^{p-1}\not\equiv 1\pmod{p}$ 必定至少占了一半，考虑证明，显然 $[1,p-1]$ 中的数与 $\bmod p$ 意义下的乘法构成了一个群，而我们如果将 $a^{p-1}\equiv 1\pmod{p}$ 的 $a$ 提取出来，那么它们与 $\bmod p$ 意义下的乘法也构成了一个群，也就是说，满足 $a^{p-1}\equiv 1\pmod{p}$ 的 $a$ 构成的群是 $[1,p-1]$ 中的数的全体构成的群的子群，其大小为 $p-1$ 的因数，又因为我们事先假定至少存在一个这样的 $a$，因此其大小至多占一半，因此对于这样的 $p$，其检验出来的概率是比较高的。

但是问题就来了，对于某些合数 $p$，不存在这样的 $a$，使得 $a^{p-1}\not\equiv 1\pmod{p}$，这样一来，我们不论随机怎样的 $a$，都不可能把这样的 $p$ hack 掉，此时用上面的方法就不太行得通了。我们考虑另一个判定质数的准则，对于质数 $p$，有满足 $x^2\equiv 1\pmod{p}$ 的 $x$ 只有 $1,p-1$，这条准则反过来说同样存在大量反例，即，存在不少合数，也有满足$x^2\equiv 1\pmod{p}$ 的 $x$ 只有 $1,p-1$。但是如果我们把上面两条结合在一起，其判定质数的概率就增加了许多，即，我们考虑这样的过程：先特判 $p$ 为偶数的情况，此时只有 $2$ 是质数，其余都是合数。对于剩余的情况那么我们假设 $p-1=2^k·t$。我们随机取 $[1,p-1]$ 中某个数 $x$ 并计算 $x^t\bmod p$，如果是 $\pm 1$ 则 hack 失败，否则我们不断平方直到平方到 $x^{(p-1)/2}$ 为止，如果某一步得到了 $p-1$ 则 hack 失败，否则则说明 $p$ 是合数。可以证明该做法单次判定质数的概率高达 $75\%$，因此你随机取 $8$ 个这样的 $x$ 出错的概率就非常小了。这就是 Miller-Rabin 素性测试。

ll _rand() {return rand() & 32767;}
ll Rand() {return (_rand() << 45) | (_rand() << 30) | (_rand() << 15) | _rand();}
typedef __int128_t i128;
const int LIM = 128;
const int PTEST[10] = {0, 2, 3, 5, 7, 13, 19, 61, 233};
inline ll mul(ll a, ll b, ll mod) {
	ll r = ((long double) a / mod * b + 0.5); r = a * b - r * mod;
	return r < 0 ? (r + mod) : r;
}
ll qpow(ll x, ll e, ll mod) {
	ll ret = 1;
	for (; e; e >>= 1, x = mul(x, x, mod))
		if (e & 1) ret = mul(ret, x, mod);
	return ret;
}
bool ptest(ll n) {
	if (n <= 500) {
		for (int i = 2; i * i <= n; i++) if (n % i == 0) return 0;
		return 1;
	} else {
		if (n % 2 == 0) return 0;
		ll tmp = n - 1;
		while (tmp % 2 == 0) tmp >>= 1;
		for (int i = 1; i <= 8; i++) {
			ll pw = qpow(PTEST[i], tmp, n); bool flg = 1;
			if (pw == 1 || pw == n - 1) flg = 0;
			while (tmp * 2 < n - 1) {
				tmp <<= 1; pw = mul(pw, pw, n);
				if (pw == n - 1) {flg = 0; break;}
			}
			if (flg) return 0;
		}
		return 1;
	}
}

现在我们已经知道了如何判定一个大数是否是质数，考虑如何对一个大数分解质因数。考虑递归，如果 $n$ 是质数则直接返回，这个可以一遍 MR 带走。否则我们考虑找到任意一个 $n$ 的真因数（即，$\ne 1$ 且 $\ne n$ 的因数）$v$，然后递归分解 $v$ 和 $\dfrac{n}{v}$，那么怎么找这个 $v$ 呢？这就要用到一个非常著名的定理，生日悖论了。根据生日悖论，期望 $23$ 人中能找到两个生日在同一天的人，更具体地，如果我们生成 $\mathcal O(\sqrt{n})$ 个 $[1,n]$ 中的数，那么期望有两个数相同。那么我们不妨假设 $n$ 存在某个因子 $D$，那么在期望条件下，如果我们生成 $\mathcal O(\sqrt{D})$ 个 $[0,n-1]$ 中的数，会存在两个数模 $D$ 同余，也就是说我们如果生成 $\mathcal O(\sqrt{D})$ 个数，那么我们期望就能够找到 $D$ 这个因子，而 $n$ 的最小质因子 $\le\sqrt{n}$，因此我们期望使用 $n^{0.25}$ 的时间就能找到 $n$ 的因子。

大体思想就是这么多，接下来就是具体实现的问题了。PR 找 $n$ 的某个真因子大致有两种实现方式：追及法和倍增法。由于 rand() 开销较大，我们不妨考虑一个随机映射 $f(x)=(x^2+C)\bmod n$，其中 $C$ 为在 PR 一开始确定的随机数，然后强制要求下一个生成的数为 $f(\text{prev})$，其中 $\text{prev}$ 为上一个生成的随机数。下面将较为详细地讲解两种实现方法：

• 追及法：任取两个数 $x_1,x_2$，然后每次令 $x_1=f(f(x_1))$，$x_2=f(x_2)$，如果某一步发现 $\gcd(n,|x_1-x_2|)\in[2,n-1]$ 就返回，根据上文，如果我们想找到因子 $D$，那么我们期望需要 $\sqrt{D}$ 步。
• 倍增法：考虑动态维护一个变量 $\text{stp}$，然后进行 $\text{stp}$ 步，每次令 $x_1=f(x_1)$，$\text{stp}$ 结束后令 $x_2=x_1$，同理，如果某一时刻 $\gcd(n,|x_1-x_2|)\in[2,n-1]$。

最后还有一点，如果直接按照上文所述实现，复杂度将会多一个 $\log$，这将导致模板题无法通过，这里有一个非常使用的 trick：设一个阈值 $\text{lim}$（一般取 $2^7=128$），然后每 $\text{lim}$ 步一个单位，每次将 $\text{lim}$ 步得到的 $|x_1-x_2|$ 乘起来，$\text{lim}$ 结束后再与 $n$ 取 $\gcd$，如果发现与 $n$ 的 $\gcd$ 大于 $1$ 再回去找具体哪一步与 $n$ 不互质。我的实现使用的是后者（也就是倍增法），即每 $128$ 步取一步 $\gcd$，$\text{stp}$ 步结束再取一个 $\gcd$，实测跑得飞快。$350$ 次 $10^{18}$ 级别的 Pollard-Rho 大约只用了 $200\text{ms}$。

下面给出模板题代码：

ll _rand() {return rand() & 32767;}
ll Rand() {return (_rand() << 45) | (_rand() << 30) | (_rand() << 15) | _rand();}
typedef __int128_t i128;
const int LIM = 128;
const int PTEST[10] = {0, 2, 3, 5, 7, 13, 19, 61, 233};
inline ll mul(ll a, ll b, ll mod) {
	ll r = ((long double) a / mod * b + 0.5); r = a * b - r * mod;
	return r < 0 ? (r + mod) : r;
}
ll qpow(ll x, ll e, ll mod) {
	ll ret = 1;
	for (; e; e >>= 1, x = mul(x, x, mod))
		if (e & 1) ret = mul(ret, x, mod);
	return ret;
}
bool ptest(ll n) {
	if (n <= 500) {
		for (int i = 2; i * i <= n; i++) if (n % i == 0) return 0;
		return 1;
	} else {
		if (n % 2 == 0) return 0;
		ll tmp = n - 1;
		while (tmp % 2 == 0) tmp >>= 1;
		for (int i = 1; i <= 8; i++) {
			ll pw = qpow(PTEST[i], tmp, n); bool flg = 1;
			if (pw == 1 || pw == n - 1) flg = 0;
			while (tmp * 2 < n - 1) {
				tmp <<= 1; pw = mul(pw, pw, n);
				if (pw == n - 1) {flg = 0; break;}
			}
			if (flg) return 0;
		}
		return 1;
	}
}
ll mx = 0;
ll _gcd(ll x, ll y) {return (!y) ? x : _gcd(y, x % y);}
ll finddiv(ll n, ll C) {
	ll c1 = Rand() % n + 1, c2 = c1;
	for (int st = 1; ; st <<= 1) {
		c2 = c1; ll prd = 1;
		for (int j = 1; j <= st; j++) {
			c1 = (mul(c1, c1, n) + C) % n;
			prd = mul(prd, abs(c1 - c2), n);
			if ((((j & 127) == 0) || j == st ) &&_gcd(prd, n) > 1)
				return _gcd(prd, n);
		}
	}
	return n;
}
void pollard_rho(ll n) {
	if (n <= mx) return;
	if (ptest(n)) {chkmax(mx, n); return;}
	ll v = finddiv(n, Rand() % (n - 1) + 1);
	pollard_rho(v); pollard_rho(n / v);
}
void solve() {
	ll n; scanf("%lld", &n);
	if (ptest(n)) {puts("Prime"); return;}
	mx = 0; pollard_rho(n); printf("%lld\n", mx);
}
int main() {
	int qu; scanf("%d", &qu);
	while (qu--) solve();
	return 0;
}

例题：

1. P5605 小A与两位神仙

首先考虑对 $m$ 分解质因数，需要 Pollard-Rho，设 $m=p^k$。

由于 $\gcd(x,m)=1,\gcd(y,m)=1$，因此我们不妨求出 $x,y$ 在模 $m$ 意义下的阶，不妨设之为 $dx,dy$，那么有结论：存在正整数 $a$ 使得 $x^a\equiv y\pmod{m}$ 当且仅当 $dy\mid dx$。考虑证明：取 $m$ 的原根 $g$，设 $x=g^{A},y=g^{B}$，那么 $dx=\dfrac{\varphi(m)}{\gcd(\varphi(m),A)},dy=\dfrac{\varphi(m)}{\gcd(\varphi(m),B)}$，那么 $x^a\equiv y\pmod{m}$ 当且仅当 $a·A\equiv B\pmod{\varphi(m)}$，根据斐蜀定理，上式有解当且仅当 $(A,\varphi(m))\mid B$，即 $(A,\varphi(m))\mid\gcd(B,\varphi(m))$，显然，这也等价于 $dy\mid dx$。

于是这题思路就非常 ez 了，直接对 $\varphi(m)$ 先分解一通质因数，然后 polylog 求阶后判一下整除性即可。

2. P6730 [WC2020] 猜数游戏

考虑转化题意：对于一组 $(i,j)$，如果存在正整数 $m$ 使得 $a_i^m\equiv a_j\pmod{p}$，则连一条 $i\to j$ 的边，那么对于一组 $S$，其答案就是最小的集合 $T$ 使得 $T\subseteq S$，且任意 $S$ 中的元素要么属于集合 $T$，要么可以从集合 $T$ 中某个元素一步到达。

由于我们连出来的这张图是一个传递闭包关系，即如果存在边 $i\to j,j\to k$，必然存在边 $i\to k$，因此，对于最优集合 $T$，必然不存在 $x,y$，使得 $x\in S,y\in T$，因此我们考虑在每个点处计算贡献，对于一个点 $x$，其会对 $S$ 集合对应的 $T$ 的大小产生 $1$ 的贡献，当且仅当其属于 $x$，但是能够到达 $x$ 的点都不属于 $S$。记能够到达 $x$ 的点数为 $cnt_x$，那么 $x$ 对答案的贡献为 $2^{n-cnt_x}$。但是这样计算实际上是错误的，因为对于同一强连通分量中的点，我们有可能需要选择其中某个点加入集合 $T$，但由于其在同一强连通分量中，所有点都是可达的，我们并不会将其计入贡献，这其实好办，我们只用为该强连通分量中的点钦定一个顺序即可。

最后考虑如何判定是否存在正整数 $m$ 使得 $a_i^m\equiv a_j\pmod{p}$，乍一看和上一题一样，但再仔细以看发现此题并没有保证 $\gcd(a_i,p)=1$，因此这题还要多一些特判——记 $p=P^k$，设 $E_i$ 表示 $a_i$ 质因数分解中 $P$ 的次数，那么如果 $E_i\ne 0,E_j\ne 0$ 则和上一题的情形一样，如果 $E_i,E_j$ 中恰有一者非零则显然不存在这样的 $m$，否则必然有 $m=\dfrac{E_j}{E_i}$，ksm check 一下即可。

时间复杂度 $\Theta(n^2)$。