学习笔记 - Pollard_Rho

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

---

Part0. 前置

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

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

---

Part1. 引入

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

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

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()，实际上采用的是伪随机数。而且是通过上一个随机化得到的数来计算出下一个数。

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

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

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。那么这个程序有没有可能死循环呢？实际上我们看到代码里这一段：

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，然后分别递归分解 m 和 num/m。

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. 源代码 <\>

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，欢迎提问。