Pollard Rho

前置芝士

1. 费马小定理

内容

对于任意质数 $p$ 和 $a$，其中 $p\nmid a$，满足：

$$a^{p-1}\equiv 1\pmod p$$

证明

根据裴蜀定理，可以知道 $\{1\dots n\}$ 与 $\{a\dots(an)\}$ 是一一对应的，因此：

$$\begin{aligned} \prod_{i=1}^{p-1}ia&\equiv\prod_{i=1}^{p-1}i&\pmod p\\ a^{p-1}\prod_{i=1}^{p-1}i&\equiv\prod_{i=1}^{p-1}i&\pmod p\\ a^{p-1}&\equiv1&\pmod p \end{aligned}$$

拓展

这个定理并没有逆定理，所以不能直接用于素性检验。但是值的注意的是，它能排除掉很大一部分合数。

2. 二次探测定理

内容

对于质数 $p$ 和 $a$，有以下性质：

如果满足：

$$a^2\equiv1\pmod p$$

则有：

$$a\equiv\pm1\pmod p$$

证明

其实很简单，就让读者自行证明了。

根据 $a^2\equiv1\pmod p$，可以知道 $a^2-1=kp$。

然后根据平方差公式可得：$kp=(a+1)(a-1)$。

因为 $p$ 是质数，所以 $a+1$ 和 $a-1$ 中，至少有一个是 $p$ 的倍数。

也就是 $a\equiv\pm1\pmod p$。

3. Miller Rabin 素性检验

算法流程

其实就是把费马小定理和二次探测定理结合了一下。

对于一个奇数 $n$，我们令 $n=2^km+1$，其中 $m$ 是奇数。

我们先随机一个数 $a$，计算出 $b=a^m$。

然后将 $b$ 平方 $k$ 次，每次平方的时候用二次探测定理判断一下。

然后经过 $k$ 次平方的 $b$ 变成了 $a^{p-1}$，这个时候就可以用费马小定理判断了。

正确率

假设你随机了 $s$ 次，那么判断错误的概率就是 $\frac{1}{4^s}$。

所以多随几次就可以了。

拓展

对于 OI 来说，我们可以考虑不随机，而是使用固定的数来素性检验。

判断 int 范围内的质数只要前 $4$ 个质数。

判断 $10^{18}$ 范围内的质数只要前 $9$ 个质数。

判断 long long 范围内的质数只要前 $10$ 个质数。

接下来给出本题范围内（即 $n\le10^{18}$）的 Miller Rabin 代码：

using ll = long long;
ll mul(ll a, ll b, ll mod){
  ll r = a * b - mod * (ll)(1.L / mod * a * b);
  return r - mod * (r >= mod) + mod * (r < 0);
}
ll qpow(ll a, ll b, ll mod) {
  ll res(1);
  for (; b; b >>= 1, a = mul(a, a, mod))
    if (b & 1) res = mul(res, a, mod);
  return res;
}
bool is_prime(ll n) {
  if (n <= 1) return false;
  vector<ll> base = {2, 3, 5, 7, 11, 13, 17, 19, 23};
  for (ll p : base) {
    if (n == p) return true;
    if (n % p == 0) return false;
  }
  ll m = (n - 1) >> __builtin_ctz(n - 1);
  for (ll p : base) {
    ll t = m, a = qpow(p, m, n);
    while (t != n - 1 && a != 1 && a != n - 1)
      a = mul(a, a, n), t *= 2;
    if (a != n - 1 && t % 2 == 0) return false;
  }
  return true;
}

值的一提的是，上述代码中的乘法没有使用 int128，而是使用了 long double，这样会极大得提升运算效率。

而且 is_prime 函数并没有完全按照之前算法流程来写，而是稍加改动，但是本质上是相同的。

4. 生日悖论

内容

在 $n$ 个数中，每次随机选取 $1$ 个数，期望 $O(\sqrt n)$ 次出现重复的数。

证明

不太会，大概用斯特林公式估算一下。

Pollard Rho

寻找一个数的真因子

首先考虑对于一个合数 $n$ 找到它的一个因子 $d$，满足 $d\mid n,1<d<n$。

大致思路

对于 $n$ 的一个奇质因子 $p$，我们把 $0$ 到 $p-1$ 中每个数 $x$ 向 $(x^2+1)\bmod p$ 连边，这样构成了一个基环森林。

不难发现这个基环森林有一些性质：

1. 首先我们可以把 $(x^2+1)\bmod p$ 看作是随机的，所以根据生日悖论，所有的环大小都是期望 $O(\sqrt p)$ 的。
2. 再观察环上挂下来的树，可以看作一棵随机的二叉树，所以期望树高是 $O(\log p)$ 的。

我们考虑两个指针 $x$ 和 $y$ ，一开始指向同一个随机节点，然后每次 $x$ 走 $1$ 步，$y$ 走 $2$ 步。

经过 $O(\log p)$ 步 $x$ 和 $y$ 都会到达环上，而根据追及问题最多再经过环的大小步，$x$ 和 $y$ 就必定会重合，也就是经过期望 $O(\sqrt p)$ 步，$x$ 和 $y$ 就会重合。

我们将 $x$ 和 $y$ 放在 $\pmod n$ 的意义下，那么就是经过期望 $O(\sqrt p)$ 步，$p\mid\gcd(n,x-y)$，此时的 $\gcd$ 还不能作为最终结果，因为还要考虑 $\gcd=n$ 的情况，但是这种情况的发生概率很小，所以直接重新选择起始点即可。

另外的，对于质因子 $2$ 可能需要特判。

这样时间复杂度就是 $O(\sqrt{\min\{p}\})$ 的，然而 $n$ 的最小质因子 $p\le \sqrt n$，所以时间复杂度也可以看作是 $O(n^{\frac14})$ 的。

一些优化

我们可以设置一个块长 $M$，把每 $M$ 步的 $x-y$ 累乘起来，一起取 $\gcd$ 判断，要注意的是，如果累乘起来的变量会变为 $0$，就不要乘上去，$M$ 一般取 $256$ 或者 $128$，但是好像没有一定要取 $2$ 的幂次这种说法。

当遇到了 $x=y$ 的情况，说明走完了模 $n$ 意义下的一大圈，其他位置都是 $\gcd(x-y,n)=1$，而只有这个位置上 $\gcd=n$，所以我们再随另一个起始点即可。

分解

直接递归，找因子前先判断一下是否是质数即可。

代码实例

/*
* @Author: ftt2333
* @Email: ftt2333@126.com
* @LastEditTime: 2022-12-21 22:18:29
*/

#include <bits/stdc++.h>
using namespace std;
#ifndef LOCAL
void debug(...) { }
#endif
#define all(v) (v).begin(), (v).end()
using ll = long long;
using lll = __int128;

mt19937_64 rnd(114514);

ll mul(ll a, ll b, ll mod){
  ll r = a * b - mod * (ll)(1.L / mod * a * b);
  return r - mod * (r >= mod) + mod * (r < 0);
}
ll add(ll a, ll b, ll mod) {
  a += b;
  return a >= mod ? a - mod : a;
}
ll sub(ll a, ll b, ll mod) {
  a -= b;
  return a < 0 ? a + mod : a;
}
ll gcd(ll a, ll b) { return b ? gcd(b, a % b) : a; }
ll qpow(ll a, ll b, ll mod) {
  ll res(1);
  for (; b; b >>= 1, a = mul(a, a, mod))
    if (b & 1) res = mul(res, a, mod);
  return res;
}
bool is_prime(ll n) {
  if (n <= 1) return false;
  vector<ll> base = {2, 3, 5, 7, 11, 13, 17, 19, 23};
  for (ll p : base) {
    if (n == p) return true;
    if (n % p == 0) return false;
  }
  ll m = (n - 1) >> __builtin_ctz(n - 1);
  for (ll p : base) {
    ll t = m, a = qpow(p, m, n);
    while (t != n - 1 && a != 1 && a != n - 1)
      a = mul(a, a, n), t *= 2;
    if (a != n - 1 && t % 2 == 0) return false;
  }
  return true;
}
ll get_factor(ll n) {
  if (n % 2 == 0) return 2;
  auto f = [&](ll x) { return add(mul(x, x, n), 1, n); };
  ll x = 0, y = 0, tot = 0, p = 1, q, g;
  for (ll i = 0; (i & 0xff) || (g = gcd(p, n)) == 1; i++, x = f(x), y = f(f(y))) {
    if (x == y) x = tot++, y = f(x);
    q = mul(p, sub(x, y, n), n);
    if (q) p = q;
  }
  return g;
}

vector<ll> solve(ll n) {
  if (n == 1) return {};
  if (is_prime(n)) return {n};
  ll d = get_factor(n);
  auto v1 = solve(d), v2 = solve(n / d);
  auto i1 = v1.begin(), i2 = v2.begin();
  vector<ll> ans;
  while (i1 != v1.end() || i2 != v2.end()) {
    if (i1 == v1.end()) ans.push_back(*i2++);
    else if (i2 == v2.end()) ans.push_back(*i1++);
    else {
      if (*i1 < *i2) ans.push_back(*i1++);
      else ans.push_back(*i2++);
    }
  }
  return ans;
}

int main() {
  cin.tie(nullptr)->sync_with_stdio(false);
  int t; cin >> t;
  while (t--) {
    ll x; cin >> x;
    auto res = solve(x);
    if (res.size() <= 1) cout << "Prime\n";
    else cout << res.back() << '\n';
  }
}