Miller Rabin

Miller Rabin 是一个快速判断质数的方法。

$a^{x-1} \equiv 1 \pmod x(1 \leqslant a \leqslant x - 1)$这个等式，在$x$是质数的时候根据费马小定理是成立的。在$x$不是质数的时候，这个等式有大概率不成立。

那么我们$\text{rand}$几个$a$，如果有至少一次不成立，就说明$x$不是质数，否则$x$有大概率是质数。多$\text{rand}$几次，正确率就会很高。

然而，对于某些合数，这个等式对大多数(甚至全部)的$a$都是成立的。于是我们需要加特技。

如果$a^b \equiv 1 \pmod x$($x$是质数且$b$是偶数)，那么$a^{\frac{b}{2}} \equiv 1 \pmod x$ 或 $a^{\frac{b}{2}} \equiv x - 1 \pmod x$。

假设$x - 1 = c \times 2^d$，如果$x$是质数：

  • 如果$a^{c \times 2^d} \equiv 1 \pmod x$，那么$a^{c \times 2^{d - 1}} \equiv 1 \pmod x$或$a^{c \times 2^{d - 1}} \equiv x - 1 \pmod x$。

  • 如果$a^{c \times 2^{d - 1}} \equiv 1 \pmod x$，那么$a^{c \times 2^{d - 2}} \equiv 1 \pmod x$或$a^{c \times 2^{d - 2}} \equiv x - 1 \pmod x$。

以此类推，如果$a^{c \times 2^{d - k}} \equiv 1 \pmod x$，那么$a^{c \times 2^{d - k - 1}} \equiv 1 \pmod x$或$a^{c \times 2^{d - k - 1}} \equiv x - 1 \pmod x$。

这些条件只要有一条不满足，那么$x$就不是质数。

据说这东西对于固定的$a$做一次的错误率约为$\frac{1}{4}$，如果多$\text{rand}$几次，正确率会很高。

Pollard Rho

Pollard Rho 是一个快速分解质因数的方法。

如果要分解$x$的质因数，先用 Miller Rabin 判断$x$是不是质数，如果是的话就不用分解了。否则我们需要找到一个数$g$使得$1 < g < x$并且$g$是$x$的约数，之后我们再递归分解$g$和$\frac{x}{g}$即可。

接下来的问题就是如何找到$g$。

先随便定义一个函数$f$，例如$f(i) = i^2 + t$($t$为任意正整数)。

你有一个正整数$p$，每次把$p$变成$f(p)$，直到$p$在模$x$意义下变成了某一个数两次后停止。那么这个操作的执行次数大约是$O(\sqrt{x})$的。

感性理解一下就是把$f$在模$x$意义下看作一个随机函数，那么大约随机$O(\sqrt{x})$次就会产生重复(和生日悖论差不多)。

如果$x$不是质数，那么一定能找到一个$g \leqslant \sqrt{x}$。进行多次把$p$变成$f(p)$的操作，那么在模$g$意义下产生重复的期望次数要明显高于在模$x$意义下重复的期望次数。

再感性理解一下，在模$g$意义下产生重复而不在模$x$意义下产生重复的期望次数约为$O(\sqrt{g}) \leqslant O (x^{\frac{1}{4}})$(因为$g \leqslant \sqrt{x}$)。

我们再考虑另一个问题，如果每次把$p$变成$f(p)$，如何快速找到在模$x$意义下发生重复的时间？(不一定要是第一次重复)

当然直接开个$set$存出现过的值是可以的，但是这里说以下另一个做法。我们开一个变量$q$，当进行操作的次数为$2$的幂次时，把$q$赋值为$p$。每次变换完，我们看一下$p$和$q$是否相等。如果$p$和$q$相等，那么我们就找到了发生重复的两个时间。

关于这个做法的复杂度，感性理解是不会超过第一次重复的时间的常数倍的，那么也是大约$O(\sqrt{x})$的。

接下来我们就讲一下 Pollard Rho 中是如何找$g$的。

一开始随机一个函数$f$和一个初始值$p$，每次操作把$p$变成$f(p)$。

再开一个变量$q$，当进行操作的次数为$2$的幂次时，把$q$赋值为$p$。

每次变换完，我们看一下$|p - q|$和$x$的最大公约数。如果这个数$ > 1$且$ < x$，那么这个数就是我们要找的$g$。

如果操作完$p$和$q$相等了，我们就找到了一次重复，那么再进行多少次操作都找不到$g$了，就重新随机$f$和$p$，从头再做一遍。

关于复杂度，如果找到了一个$g$，那么说明它在模$g$意义下产生了重复并且没有在模$x$意义下产生重复，根据上面的分析(感性理解)，是$O(x ^ {\frac{1}{4}})$的。

由于在分解质因数过程中，上面的过程会被调用$O(\log{x})$次，所以总复杂度是$O(x^{\frac{1}{4}}\log{x})$(忽略快(慢)速乘的复杂度)。

题目链接(垃圾题目卡我常数，最后对着数据打表才过)

以下代码删去了输入输出和打表：

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
75
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;

namespace millerrabin {

    const int testcnt = 7;

    inline ll add(ll x, ll y, ll mod) {
        x += y; return x < mod ? x : x - mod;
    }

    ll mul(ll x, ll y, ll mod) {
        ll ans = 0;
        while (y) {
            if (y & 1) ans = add(ans, x, mod);
            y >>= 1; x = add(x, x, mod);
        }
        return ans;
    }

    ll fpow(ll x, ll y, ll mod) {
        ll ans = 1;
        while (y) {
            if (y & 1) ans = mul(ans, x, mod);
            y >>= 1; x = mul(x, x, mod);
        }
        return ans;
    }

    bool main(ll x) {
        if (x == 1) return 0;
        for (int i = 1; i <= testcnt; ++i) {
            ll v = rand() % (x - 1) + 1, now = fpow(v, x - 1, x);
            if (now != 1) return 0;
            for (ll t = x - 1; ~t & 1; t >>= 1) {
                now = fpow(v, t >> 1, x);
                if (now == x - 1) break;
                else if (now != 1) return 0;
            }
        }
        return 1;
    }
}

namespace pollardrho {
    using millerrabin::add;
    using millerrabin::mul;
    vector<ll> ans;

    void solve(ll x) {
        if (millerrabin::main(x)) ans.push_back(x);
        else {
            while (1) {
                ll p = rand() % (x - 1) + 1, q = p, d = rand() % (x - 1) + 1;
                for (int i = 1; ; ++i) {
                    p = add(mul(p, p, x), d, x);
                    if (p == q) break;
                    ll g = __gcd(abs(p - q), x);
                    if (g > 1 && g < x) {
                        solve(g); solve(x / g);
                        return;
                    }
                    if ((i & -i) == i) q = p;
                }
            }
        }
    }

    vector<ll> main(ll x) {
        ans.clear();
        if (x > 1) solve(x);
        return ans;
    }
}