「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);
}
}