Pollard-Rho 算法可以用期望 $O(\sqrt[4]{n})$ 的时间找到合数 $n$ 的一个非平凡因子。

算法流程

这里讲解改进后的 Pollard-Rho 算法变种。

在一次 Pollard-Rho 中，我们随机一个数 $c$，并设定阈值 $k$，初始为 $2$，每次循环后加倍。

对于每一次循环，令 $y$ 等于当前的 $x$，然后进行 $k$ 次，每次让 $x \leftarrow f(x) = (x ^ 2 + c) \bmod n$，然后检查 $\vert x - y \vert$ 和 $n$ 是否有非平凡的 $\gcd$，如果有则退出循环。

实践后发现这样速度较快，这可说明该算法的复杂度较低，但是笔者不会严格证明算法的复杂度。有一个小优化，就是每次不直接检查 $\gcd$，而是将变量 $z \leftarrow z \times \vert x - y \vert$。进行 $B$ 次乘法后只要检查一次 $z$ 和 $n$ 的 $\gcd$ 即可。这可以节省 $\gcd$ 使用的时间，$B$ 一般使用 $128$。

Pollard-Rho 结合 Miller-Rabin 可以实现快速素因数分解。

代码实现

「模板」Pollard-Rho 算法（Luogu 4718）

1	#include <ctime>
2	#include <cstdio>
3	#include <cstdlib>
4
5	typedef long long ll;
6	int T;
7	ll n, ans;
8
9	inline ll mult(const ll &a, const ll &b, const ll &m) {
10	return (__int128) a * b % m;
11	}
12
13	inline ll func(const ll &x, const ll &n) {
14	return x < n ? x : x - n;
15	}
16
17	ll qpow(ll a, ll b, ll m) {
18	ll c = 1;
19	for (; b; b >>= 1, a = mult(a, a, m)) {
20	if (b & 1) c = mult(a, c, m);
21	}
22	return c;
23	}
24
25	bool miller(ll m, ll d, ll r, ll n) {
26	if (m > n - 2) return true;
27	ll x = qpow(m, d, n);
28	if (x == 1 \|\| x == n - 1) return true;
29	for (int i = 0; i < r; i++) {
30	x = mult(x, x, n);
31	if (x == n - 1) return true;
32	}
33	return false;
34	}
35
36	bool is_prime(ll n) {
37	if (n <= 2) return n == 2;
38	static ll m[] = { 2, 325, 9375, 28178, 450775, 9780504, 1795265022 };
39	ll d = n - 1, r = 0;
40	while (~d & 1) d /= 2, r++;
41	for (int i = 0; i < 7; i++) {
42	if (!miller(m[i], d, r, n)) return false;
43	}
44	return true;
45	}
46
47	ll randl() {
48	return rand() \| ((ll) rand() << 31);
49	}
50
51	ll gcd(ll a, ll b) {
52	return b ? gcd(b, a % b) : a;
53	}
54
55	ll absl(ll x) {
56	return x < 0 ? -x : x;
57	}
58
59	ll factor(ll n) {
60	ll c = randl() % (n - 1) + 1, x = 0, y = 0, d = 1;
61	for (ll k = 2, t = 1; ; k *= 2, y = x, t = 1) {
62	for (ll i = 1; i <= k; i++) {
63	x = func(mult(x, x, n) + c, n), t = mult(t, absl(x - y), n);
64	if (!(i & 127) && (d = gcd(t, n)) > 1) break;
65	}
66	if (d > 1 \|\| (d = gcd(t, n)) > 1) break;
67	}
68	return d;
69	}
70
71	void solve(ll n) {
72	if (n == 1 \|\| n <= ans) return;
73	if (is_prime(n)) {
74	ans = n;
75	return;
76	}
77	ll x = factor(n);
78	while (x == n) x = factor(n);
79	while (n % x == 0) n /= x;
80	solve(x), solve(n);
81	}
82
83	int main() {
84	scanf("%d", &T);
85	srand(time(0) ^ T);
86	while (T--) {
87	scanf("%lld", &n);
88	ans = 0;
89	solve(n);
90	if (ans == n) {
91	puts("Prime");
92	} else {
93	printf("%lld\n", ans);
94	}
95	}
96	return 0;
97	}