Pollard-Rho算法
/ / 阅读耗时 12 分钟 本文谈一个数论方面的玄学算法。
Miller-Rabin素数测试
平时判定一个数是不是质数，可以采用试除法，在$O(\sqrt n)$下解决。另外可以在$O(n)$线性筛出质数，这两种方法对空间和时间有巨大的消耗。
有一种方法可以很快地判定一个数是不是质数，这就是MillerRabin素数测试法，它虽然有误判的概率，不过在平时使用中，准确度可以达到要求。
首先需要知道一个定理：费马小定理：
【费马小定理】若$p$是质数，且$a$与$p$互质，那么有$a^{p-1}\equiv 1\pmod p$
它给出了一个质数的必要条件，不过不是充分条件。当然，如果对于一个给定的$a$，$a^{p-1}$并不与1同余，那么这个$p$一定不是质数，这就是费马测试。
费马测试当然不能充分地判定质数，有许多合数同样可以通过费马测试，因此还需要另一个东西来判定。这里用二次探测定理：
【二次探测定理】若$p$是质数，$a^2\equiv 1\pmod p$那么$a\equiv 1\pmod p$或者$a\equiv p-1\pmod p$
这里可以与费马小定理同时使用来判定质数。
Miller-Rabin测试的步骤是这样的。首先选定若干个质数，比如这里选择2、3、7、61、24251这五个，然后用这五个数依次进行判定。对于选定的一个数$a$，求$a^{p-1}$对$p$的模，根据费马小定理初步判定是否为素数，若不为1，则肯定为合数，判定完毕。
如果为1，那么检测$p-1$是否为偶数，若是，那么$p-1$可以写成$p-1=2^km$的形式，那么就可以写成：
就是：
然后将$p-1$除以2，用二次探测定理向下依次判别，直到不符合二次探测定理或者p-1成为奇数或者不符合二次探测定理的使用条件。
经过这五轮测试后，判定质数的正确性是相当高的。在$10^{16}$的数据中，只有一个数$46856248255981$作为一个合数却通过了所有测试。
下面给出代码，判定给定的n个数哪些是质数：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
using namespace std;
const int P[] = {2, 3, 7, 61, 24251};
long long qPow(long long a, long long b, long long mod) {
__int128 ans = 1, sta = a;
while (b) {
if (b & 1)ans = ans * sta % mod;
sta = sta * sta % mod, b >>= 1;
}
return static_cast<long long int>(ans);
}
bool mr(long long x, int s) {
if (qPow(P[s], x - 1, x) != 1)return false;
long long k = x - 1;
while ((k & 1) == 0) {
k >>= 1;
long long d = qPow(P[s], k, x);
if (d != 1 && d != x - 1)return false;
if (d == x - 1)return true;
}
return true;
}
bool MR(long long x) {
if (x == 2 || x == 3 || x == 7 || x == 61 || x == 24251)return true;
return mr(x, 0) & mr(x, 1) & mr(x, 2) & mr(x, 3) & mr(x, 4);
}
int main() {
int n;
scanf("%d", &n);
long long s;
while (scanf("%lld", &s) != EOF) {
if (MR(s))printf("Prime\n");
else printf("Not prime\n");
}
return 0;
}
Pollard-Rho算法
现在我们能很快地判定质数了，但是如果想分解质因数，是否有更快的解法？
普通的分解质因数方法也是试除，复杂度在$O(\sqrt n)$，PR算法能将它优化到$O(n^{0.25})$（？），其实是玄学复杂度。
思想比较简单粗暴：在$1~n$中选很多随机数，然后检测它们是不是$n$的因子。
这样的话概率可能太低，根本无法做到快速判别。于是Pollard提出了一种生成随机数的方法。首先先选定一个随机数$c$和一个起始值$v_0$，然后我们用这个式子迭代求出很多数：
然后用相邻两个数差的绝对值作为待检测的因子。
这样会有一个问题：这个迭代过程最终会出现循环！这是由于每一个数的生成完全依赖于前一个，并且在模意义下生成，一定会出现循环。下面会介绍一种方法解决这个问题。
对于相邻的两个数的绝对值，难道直接去试试$n$能否整除这个数？这样做概率还是太低了，可以考虑另一种思路，就是求这个数与$n$的最大公因数（$gcd$），如果最大公因数不为1，那么这个最大公因数就是$n$的一个因子。
每产生一个数就要求一次$gcd$，效率太低，不妨先迭代出很多次，将这些差的绝对值乘到一起，并模$n$，根据数论上的内容，这样做不会影响其与$n$的最大公因数，因此是可行的。
于是迭代一段时间后同一进行一次$gcd$运算，效率有了提升。不过多长时间后才进行一轮运算？可以选择倍增的方式，第一次迭代1次，之后按2、4、8的形式去扩大范围，当在某一段中发现$v$与初始的$v_0$相等，则退出（每一轮都要更新$v_0$），表示分解失败，这时候重新进行一次即可。
质数没法分解质因数，故在此之前先用Miller-Rabin测试先判一下。
下面给出洛谷模板题的示例代码：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
76
77
78
79
80
81
using namespace std;
const int P[] = {2, 3, 7, 61, 24251};
typedef long long ll;
ll qPow(ll a, ll b, ll mod) {
__int128 ans = 1, sta = a;
while (b) {
if (b & 1)ans = ans * sta % mod;
sta = sta * sta % mod, b >>= 1;
}
return ans;
}
bool mr(ll x, int s) {
if (qPow(P[s], x - 1, x) != 1)return false;
ll k = x - 1;
while ((k & 1) == 0) {
k >>= 1;
ll d = qPow(P[s], k, x);
if (d != 1 && d != x - 1)return false;
if (d == x - 1)return true;
}
return true;
}
bool MR(ll x) {
if (x == 46856248255981ll)return false;
if (x == 2 || x == 3 || x == 7 || x == 61 || x == 24251)return true;
return mr(x, 0) & mr(x, 1) & mr(x, 2) & mr(x, 3) & mr(x, 4);
}
ll gcd(ll a, ll b) {
return b == 0 ? a : gcd(b, a % b);
}
inline ll f(__int128 x, __int128 c, __int128 n) {
return (x * x + c) % n;
}
ll PR(ll x) {
ll c = rand() % (x - 1) + 1, v0 = rand() % x, k = 1, pk = 0;//生成随机数
__int128 tmp = 1;
ll v = v0;
while (true) {
v = f(v, c, x), tmp = tmp * abs(v - v0) % x;
if (tmp == 0 || v == v0)return x;
if (++pk == k) {
ll d = gcd(tmp, x);
if (d > 1)return d;
v0 = v, k <<= 1;
}
}
}
ll getAns(ll s) {//获得最大质因子
if (s == 1 || MR(s))return s;//为1或者是质数，直接返回
ll ss = s;
while (ss == s)ss = PR(s);//找出一个因子
while (s % ss == 0)s /= ss;
return max(getAns(ss), getAns(s));
}
int main() {
// freopen("text.in", "r", stdin);
srand(time(0));
int t;
cin >> t;
while (t--) {
long long s;
scanf("%lld", &s);
if (MR(s))printf("Prime\n");
else printf("%lld\n", getAns(s));
}
return 0;
}