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)\)。
浙公网安备 33010602011771号