「Pollard-Rho」学习笔记

发表于 | 分类于 |

Pollard-Rho 是用来分解大数质因子的一个玩意。

Miller-Rabin 素性判断

众所周知，费马告诉我们对于一个素数 \(p\)，有 \(a^{p - 1} \equiv 1 \pmod p\)。
然鹅它的逆命题是不成立的，即 \(a^{p - 1} \equiv 1 \pmod p\) 不是 \(p\) 的充分条件。

不过，大部分合数都可以被这个式子筛掉，所以我们可以基于这个式子进行下一步改进。
然后有二次探测定理：对于素数 \(p\)，如果有 \(a^2 \equiv 1 \pmod p\)，那么有 \(a \equiv 1 \pmod p\) 或 \(a \equiv -1 \pmod p\)。

简单地证明一下：
\(x^2 = (x + 1)(x - 1) + 1\)。
所以 \((x + 1)(x - 1) \equiv 0 \pmod p\)。
要么 \(x + 1 \equiv 0 \pmod p\)，要么 \(x - 1 \equiv 0 \pmod p\)。
然后易得那两个式子。

考虑到费马小定理和二次探测定理的形式很像，并且 \(p - 1\) 一定为偶数，所以我们可以一直运用二次探测定理。
对于不同的 \(a\)，我们把上述算法称之为 以 \(a\) 为底的 Miller-Rabin 测试。

一般底数仍然是随机选取，但当待测数不太大时，选择测试底数就有一些技巧了。
比如，如果被测数小于 \(4759123141\)，那么只需要测试三个底数 \(2,7,61\) 就足够了。
当然，你测试的越多，正确的范围肯定也越大。如果你每次都用前 \(7\) 个素数 \((2, 3, 5, 7, 11, 13, 17)\) 进行测试，所有不超过 \(341550071728320\) 的数都是正确的。
如果选用 \(2, 3, 7, 61, 24251\) 作为底数，那么 \(10^{16}\) 内唯一的强伪素数为 \(46856248255981\)。
这样的一些结论使得 Miller-Rabin 算法在OI中非常实用。
——Matrix67

根据本人测试，对于上面那个 \(46856248255981\)，增加底数 \(19****17\) 即可。
（大雾）

代码

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
//洛谷 3383.「模板」线性筛素数
#include <cstdio>
using namespace std;
int n,m;
inline int fpow(int a,int b,int mod)
{
    int ret = 1;
    for(;b;b >>= 1)
        (b & 1) && (ret = (long long)ret * a % mod),a = (long long)a * a % mod;
    return ret;
}
inline bool witness(int x,int p)
{
    if(fpow(x,p - 1,p) != 1)
        return 0;
    int k = p - 1;
    while(!(k & 1))
    {
        k >>= 1;
        int t = fpow(x,k,p);
        if((t ^ 1) && (t ^ p - 1))
            return 0;
        if(t == p - 1)
            return 1;
    }
    return 1;
}
inline bool Miller_Rabin(int n)
{
    const int base[] = {2,3,7,61,24251,19260817};
    if(n == 1)
        return 0;
    for(register int i = 0;i < 6;++i)
        if(n == base[i])
            return 1;
        else if(!witness(base[i],n) || !(n % base[i]))
            return 0;
    return 1;
}
int main()
{
    scanf("%d%d",&n,&m);
    while(m--)
    {
        scanf("%d",&n);
        puts(Miller_Rabin(n) ? "Yes" : "No");
    }
}

Pollard-Rho

考虑对于一个合数 \(n\)，它的一个质因子为 \(p\)，如果有 \(a \equiv b \pmod p\) 但 \(a \not \equiv b \pmod n\)，那么 \(|a - b|\) 便是 \(n\) 的一个非平凡因子。
所以一个比较暴力的方法就是一直随机两个数并计算它们的差与 \(n\) 的 \(\gcd\) 然后递归做。

然鹅这样子考虑到随机数的效率和质量，并没有什么卵用……
我们考虑手写随机数，设我们生成的第 \(i\) 个随机数是 \(x_i\)，并选取一个常数 \(c\)，有 \(x_i \equiv x_{i - 1}^2 + c \pmod n\)。
但是接下来又有一个问题，怎样提高准确率？
显然地，\(x\) 数列会出现循环。如果 \(x_i \equiv x_j \pmod d\ (i < j)\)，那么数列在模 \(d\) 意义下是周期的，其长度整除 \(j - i\)。
那么对于 \(j - i\) 最小的 \(> i\) 的倍数 \(s\)，有 \(x_s \equiv x_{2s} \pmod d\)。
所以我们可以对于 \(x_k\) 与 \(x_{2k}\) 进行这个过程。

又有一个问题，如果判环？
哈希？太麻烦了！
考虑在一个环上，甲、乙从同一个点出发，甲经过一条边的时间乙能经过两条边，甲、乙向同一个方向走。
那么当乙追上甲时，甲一定至少走了一圈。
这个方式恰好可以与上面的过程结合起来。

在实际应用中，一定要注意先用线性筛出空间允许的范围内的素数并试除，最后再上 Pollard-Rho。
因为它只能找到部分较大的质因子。

代码

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
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
//洛谷 4718.「模板」Pollard-Rho 算法
#include <cstdio>
#include <cmath>
#include <algorithm>
using namespace std;
long long ans;
inline long long fmul(long long a,long long b,long long mod)
{
    return (a * b - (long long)((long double)a / mod * b) * mod + mod) % mod;
}
inline long long fpow(long long a,long long b,long long mod)
{
    long long ret = 1;
    for(;b;b >>= 1)
        (b & 1) && (ret = fmul(ret,a,mod)),a = fmul(a,a,mod);
    return ret;
}
inline bool witness(long long x,long long p)
{
    if(fpow(x,p - 1,p) ^ 1)
        return 0;
    long long k = p - 1;
    while(!(k & 1))
    {
        k >>= 1;
        long long t = fpow(x,k,p);
        if((t ^ 1) && (t ^ p - 1))
            return 0;
        if(t == p - 1)
            return 1;
    }
    return 1;
}
inline bool Miller_Rabin(long long n)
{
    if(n == 1)
        return 0;
    const int base[] = {2,3,7,61,24251,19260817};
    for(register int i = 0;i < 6;++i)
        if(n == base[i])
            return 1;
        else if(!witness(base[i],n))
            return 0;
    return 1;
}
void Pollard_Rho(long long n)
{
    if(n <= ans)
        return ;
    if(Miller_Rabin(n))
    {
        ans = max(ans,n);
        return ;
    }
    const long long c = 20060312;
    const long long k = 127;
    long long t1 = 20070119;
    long long t2 = (fmul(t1,t1,n) + c) % n;
    long long i = 0,p = 1;
    while(t1 ^ t2)
    {
        ++i,p = fmul(p,abs(t1 - t2),n);
        if(!p)
        {
            long long d = __gcd(abs(t1 - t2),n);
            if((d ^ 1) && (d ^ n))
            {
                Pollard_Rho(d),Pollard_Rho(n / d);
                return ;
            }
        }
        if(!(i % k))
        {
            long long d = __gcd(p,n);
            p = 1;
            if((d ^ 1) && (d ^ n))
            {
                Pollard_Rho(d),Pollard_Rho(n / d);
                return ;
            }
        }
        t1 = (fmul(t1,t1,n) + c) % n,t2 = (fmul(t2,t2,n) + c) % n,t2 = (fmul(t2,t2,n) + c) % n;
    }
    long long d = __gcd(p,n);
    if((d ^ 1) && (d ^ n))
        Pollard_Rho(d),Pollard_Rho(n / d);
}
int T;
long long n;
int main()
{
    scanf("%d",&T);
    while(T--)
    {
        ans = 0;
        scanf("%lld",&n);
        Pollard_Rho(n);
        ans == n ? puts("Prime") : printf("%lld\n",ans);
    }
}