学习笔记 - Pollard_Rho

这是第 n 次恶补数论了……
预估在我 AFO 之前还会恶补数论 $∞$ 次 QwQ


Part0. 前置

你还得先学会这些……然后开始递归学习

Miller_Rabin 素性测试,这个我自己也有写 blog,reader 们可以去看一下 awa。


Part1. 引入

我们常会处理这样一个问题——将 $n$ 分解质因数。毕竟分解质因数过后就可以玩一些类似于组合数、贪心之类的东西,与欧拉函数 $\phi$ 也或多或少有些关系。

如果 $n$ 不是很大,意思是可以接受 $O(\sqrt n)$ 的复杂度。我们有一个非常暴力的解决方法 —— 如果 $n$ 本身不是一个质数,就一定可以把它分解成一个小于等于 $\sqrt n$ 的质数与另一个数的乘积。因此我们可以这样分解:

1
2
3
4
5
6
7
8
int function(int n){
for(int i=2;i*i<=n;i++)
while(n%i==0){
n/=i;
//i是一个质因数
}
if(n>1) /*剩下的n是一个质因数*/;
}

显然这是一个 $O(\sqrt n)$ 的算法。

但是有一些题非常恶心心……把 $n$ 的大小开到 $10^{18}$ 的范围,或者是让你分解很多次质因数。这样的算法复杂度无法解决。

于是就出现了 Pollard_Rho 算法。


Part2. 初步&引申

实际上 Pollard_Rho 的算法思想是基于这样一个发现:

【问题】

在 $1$ ~ $n$ 的数内随机地找一个元素 $a$。

如果我们直接 rand() 来找这个 $a$ ,期望复杂度为 $O(n)$。效率非常之低……

但是如果每次随机化两个数 x,y ,然后计算 x-y,因为 x-y=a 的数量要比直接随机化 x=a 的数量要多,我们可以证明它的期望复杂度是 $O(\sqrt n)$。

(详细论证见 生日悖论

再看我们现在的问题——要找一个数 $x$ 使得 $x$ 是 $n$ 的一个因数(还不是求质因数)。

最简单的随机化的思路就是在 $[1,n]$ 中随机化一个 $x$,使得 $x|n$ 。
然后根据上面提到的思路,我们可以每次随机化两个数 $x,y$,使得 $|x-y|$ 是 $n$ 的因数。

但是这样的效率仍然太低。其实不难发现,我们可以把问题这样等价转换——找一个数 $x$ 使得 $gcd(x,n)\not= 1$。

于是我们发现可以进一步提高效率:每次随机化两个数 $x,y$,使得 $gcd(|x-y|,n)\not=1$。那么我们求得的不为 $1$ 的 $gcd(|x-y|,n)$ 就是 $n$ 的一个因数。

比较显然的,满足“$gcd(x,n)\not=1$ 的 $x$” 比 满足“$x|n$ 的 $x$” 多得多(至少不会比第二种少)。

所以我们找到这样的 $x$ 的概率就会高一些。


Part3. Pollard_Rho

虽然说着是随机化 rand(),实际上采用的是伪随机数。而且是通过上一个随机化得到的数来计算出下一个数。

1
2
//x 是当前得到的一个随机数,现在要根据 x 计算出下一个“随机”数
x = ( x*x+c ) % m;

于是 Pollard_Rho 的主程序就出现了:

1
2
3
4
5
6
7
8
9
10
11
long long Pollard_Rho(long long num){
long long A = rand()%(num-2)+2 ; //在 [2,num-1] 之间取一个数
long long B = A , C = rand()%(num-1)+1 , Vgcd = 1 ;
while(Vgcd==1){ //直到不为1才退出
A = ( A*A+C ) % num;
B = ( B*B+C ) % num;
B = ( B*B+C ) % num; //B作两次,A作一次,那么 AB 就可以相对随机
Vgcd = GCD( ABS(A-B),num );
}
return Vgcd;
}

显然这样随机化一定会产生一个环,即从 $x$ 开始,进行多次 $x=(x^2+c)\bmod m$ 后一定会回到 $x$ 。

如果运气不好,可能绕一圈过后,都没能找到一个不为 1 的 gcd。那么这个程序有没有可能死循环呢?实际上我们看到代码里这一段:

1
2
3
A = ( A*A+C ) % num;
B = ( B*B+C ) % num;
B = ( B*B+C ) % num;

其实是非常重要的,如果 GCD(ABS(A-B),num) 一直为 1,经过一次循环后,A 和 B 会相等,然后 A-B 会等于 0。众所周知 GCD(0,num)=num ,所以这个时候 while 循环就会退出来。并且 Pollard_Rho 返回 num(表示查找失败)。


Part4. 分解质因数

其实也看到了,Pollard_Rho 的主函数只是分解出一个 num 的因数而已,我们要把 num 分解质因数!

这里可以想到用递归。

如果 num 是质数(这里就要用 Miller_Rabin 了),那么它本身就是一个质因数,否则用 Pollard_Rho 找到一个因数 m,然后分别递归分解 mnum/m

1
2
3
4
5
6
7
8
9
void DivideNum(long long num){
if(num==1) return; //别忘了
if(CheckPrime(num)){ //这个就用 Miller_Rabin 了
//找到一个质因数
return;
}
long long m=Pollard_Rho(num);
DivideNum(m);DivideNum(num/m);
}

Part5. 源代码 <\>

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
typedef long long ll;
typedef long double ld;
typedef unsigned long long ull;
const ll Cc[]={2,3,5,7,11,13,17,19,23,29};

ll Mul(ll Va,ll Vb,ll Vm){
ll Vc=(ld)Va/Vm*Vb;
ll Rr=(ull)Va*Vb-(ull)Vc*Vm;
return (Rr+Vm)%Vm;
}
ll QPow(ll Va,ll Vb,ll Vm){
ll Rr=1;
while(Vb){
if(Vb&1) Rr=Mul(Rr,Va,Vm);
Va=Mul(Va,Va,Vm);
Vb>>=1;
}
return Rr;
}
bool MillerRabin(ll Va,ll Vd,ll Vn){
ll Vx=QPow(Va,Vd,Vn);
if(Vx==1 || Vx==Vn-1) return true;
while(Vd!=Vn-1){
Vd*=2;Vx=Mul(Vx,Vx,Vn);
if(Vx==1) return false;
if(Vx==Vn-1) return true;
}
return false;
}
bool CKPr(ll Vn){
for(int i=0;i<10;i++){
if(Vn==Cc[i]) return true;
if(Vn%Cc[i]==0) return false;
}
ll Vd=Vn-1;
while(Vd%2==0) Vd/=2;
for(int i=0;i<10;i++)
if(!MillerRabin(Cc[i],Vd,Vn))
return false;
return true;
}
ll GCD(ll Va,ll Vb){return Vb? GCD(Vb,Va%Vb):Va;}
ll ABS(ll Vx){return Vx<0? -Vx:Vx;}
ll PollardRho(ll Vn){
for(int i=0;i<10;i++)
if(Vn%Cc[i]==0)
return Cc[i];
ll Va=1ll*rand()*rand()%(Vn-2)+2,Vb=Va,
Vg=1,Vc=1ll*rand()*rand()%(Vn-1)+1;
while(Vg==1){
Va=(Mul(Va,Va,Vn)+Vc)%Vn;
Vb=(Mul(Vb,Vb,Vn)+Vc)%Vn;
Vb=(Mul(Vb,Vb,Vn)+Vc)%Vn;
Vg=GCD(ABS(Va-Vb),Vn);
}
return Vg;
}
ll MinP(ll Vn){ //这道题是要求最小的质因数
if(Vn==1) return (1ll<<60);
if(CKPr(Vn)) return Vn;
ll Vm=PollardRho(Vn);
return min(MinP(Vm),MinP(Vn/Vm));
}

The End

Thanks for reading!

Email: lucky_glass@foxmail.com,欢迎提问。

0%