简介

Pollard Rho 算法利用生日悖论，可以在优秀的复杂度内找到一个合数的非 11 非自身因数，使用时需要搭配 Miller Rabin 算法首先判断一个数是否为质数。

生日悖论

关于生日悖论，就是说 kk 个人中两个人生日相同的概率实际上很高，一年有 nn 天，这个概率是：

P=1−∏i=0k−1n−in=1−n!nk(n−k)!P=1-\prod_{i=0}^{k-1}\frac{n-i}{n}=1-\frac{n!}{n^k(n-k)!}

当 k=22,n=365k=22,n=365 时，已经有 P≈0.476P\approx 0.476，举这个例子只是说明它的增长很快。

可以证明，随机地从 1∼n1\sim n 中选择数字，O(n)O(\sqrt n) 次会选到重复的数字。所以，当一个数 nn 是合数，它的最大质因子 pp 是 O(n)O(\sqrt n) 级别的，在模 pp 的意义下找到两个相同的数字的复杂度是 O(p)O(\sqrt p) 即 O(n14)O(n^{\frac{1}{4}})。

原理

设 nn 的一个因子为 p∈(1,n)p\in(1,n)，随机取 c∈[1,n)c\in[1,n)，定义序列 xi+1=(xi2+c) mod nx_{i+1}=(x_i^2+c)\bmod n 并令 x1=0x_1=0，可以得到一个比较随机的 0∼n−10\sim n-1 的序列，根据鸽巢原理 xix_i 是一定存在环的，并且在模 pp 的意义下环的长度是 O(p)O(\sqrt p)。

若 p∣xj−xip\mid x_j-x_i，那么就有 (xj−xi,n)=kp(x_j-x_i,n)=kp，并且 xj+1−xi+1≡(xj−xi)(xj+xi)(modn)x_{j+1}-x_{i+1}\equiv (x_j-x_i)(x_j+x_i)\pmod n，也就是说 p∣xj+1−xi+1p\mid x_{j+1}-x_{i+1}，于是只要改变 i,ji,j 的相对大小就可以判断出环上间距为 j−ij-i 的 gcd⁡\gcd 都是 pp 的倍数 。

期望找到环的复杂度是 O(p)O(\sqrt p)，并且同时用快慢指针，也是在 O(p)O(\sqrt p) 的复杂度内找到 (i,j)(i,j) 使得 p∣xj−xip\mid x_j-x_i，加上求 gcd⁡\gcd 的复杂度，算法的期望复杂度是 O(n14log⁡n)O(n^{\frac{1}{4}}\log n)。

实现

输入不超过 350350 个数字，对于每个数字 2≤n≤10182\le n\le 10^{18} 检验是否是质数，是质数就输出 Prime；如果不是质数，输出它最大的质因子是哪个。

题目链接：P4718。

先用 Miller Rabin 判掉素数，然后再用 Pollard Rho 找因数，复杂度 O(n14log⁡2n)O(n^{\frac{1}{4}}\log^2 n)，第二个 log⁡\log 代表质因数的数量。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
#include <bits/stdc++.h>
using namespace std;

namespace IO {}
using IO::read;
using IO::write;

#define int long long
mt19937_64 rng(chrono::steady_clock::now().time_since_epoch().count());

int gcd(int a, int b) {
    if (!b) return a;
    return gcd(b, a % b);
}

int qmi(int a, int k, int p) {
    int res = 1;
    while (k) {
        if (k & 1) res = (__int128)res * a % p;
        a = (__int128)a * a % p;
        k >>= 1;
    }
    return res;
}

int rho(int n) {
    if (n % 2 == 0) return 2;
    int c = uniform_int_distribution<int>(1, n-1)(rng);
    auto f = [=](int x){ return ((__int128)x * x + c) % n; };
    int x = 0, y = f(x);
    while (x != y) {
        int d = gcd(abs(x - y), n);
        if (d > 1) return d;
        x = f(x), y = f(f(y));
    }
    return n;
}

bool test(int a, int u, int k, int n) {
    int v = qmi(a, u, n);
    if (v == 1 || v == 0 || v == n-1) return true;
    for (int i = 1; i <= k; i++) {
        v = (__int128)v * v % n;
        if (v == n-1 && i != k) return true;
        if (v == 1) return false;
    }
    return false;
}

bool mr(int n) {
    static int bases[] = {2, 325, 9375, 28178, 450775, 9780504, 1795265022};  
    if (n < 3 || n % 2 == 0) return n == 2;
    int u = n-1, k = 0;
    while (u % 2 == 0) u /= 2, k++;  
    for (int a: bases) if (!test(a, u, k, n)) return false;
    return true;
}

int solve(int n) {
    static map<int, int> mp;
    if (mp.count(n)) return mp[n];

    if (mr(n)) return mp[n] = n;
    int x = n;
    while (x == n) x = rho(n);
    return mp[n] = max(solve(x), solve(n / x));
}

signed main() {
    int n, x, t;
    read(t);
    while (t--) read(n), ((x = solve(n)) == n) ? write("Prime") : write(x);
    return 0;
}