关于Miller-Rabin与Pollard-Rho算法的理解（素性测试与质因数分解）

前置

费马小定理（即若P为质数，则(A^Pequiv A pmod{P})）。
欧几里得算法（GCD）。
快速幂，龟速乘。

素性测试

引入

素性测试是OI中一个十分重要的事，在数学毒瘤题中有着举足轻重的地位。
常见的素性测试如下：

int check(int N){
	for(int i=2;i*i<=N;i++)
		if(N%i==0)return 0;
	return 1;
}

以上是一个(O(sqrt{N}))的算法，虽然不优，但在绝大多数情况下是可以的。
但是，假如(N)的范围达到了(1e18)，以上算法很明显就不行了，我们得考虑更优的算法。
引入Miller-Rabin算法。

定理1

定理1：如果(P)为一个大于2的素数，那么方程(X^2equiv1 pmod{P})的解只有1或者-1。

（是真得水的） 证明如下：
由(X^2equiv1 pmod{P})，得((X^2-1)equiv0 pmod{P})，
即({(X+1)(X-1)}equiv0 pmod{P})。
则(P mid (X+1))或(Pmid(X-1))，
即((X+1)equiv0 pmod{P})或((X-1)equiv0 pmod{P})。
则(Xequiv-1 pmod{P})或(Xequiv1 pmod{P})，即证。
（注:由于P>2，不能同时满足(Xequiv-1 pmod{P})与(Xequiv1 pmod{P})。）

定理1的延伸

定理1的逆否定理即为：
如果有(P>2)且(X^2equiv1 pmod{P})且(X otin{1,-1})，则(P)不为质数。

证明如下：
假设(P)为质数，则(P)满足定理1，即方程(X^2equiv1 pmod{P})的解只有1或者-1。
与条件矛盾，故(P)不为质数。

Miller-Rabin

基于以上定理，我们可以开始了。
假如我们要验证(P(P>2))是否是一个质数。

先设(P)是一个素数。
那么对于所有满足(A^2equiv1 pmod{P})的整数(A)，都有(Ain{1,-1})。
但如果我们枚举(A)来暴力判断，则复杂度又上升到了(O(N))。
（不妨令(A<P)，以下(A)均小于(P)）

考虑做判断时，如何快速的找到满足(A^2equiv1 pmod{P})的整数(A)。
我们考虑将((P-1))分解成(2^S·D)的形式，其中(D)为奇数。
（由于(P>2)且(P)为素数，则((P-1))为偶数，即(Sge1)）
则根据费马小定理，有

[A^{P-1}equiv 1pmod{P} ]

[A^{2^S·D}equiv 1pmod{P} ]

[({A^{2^{S-1}·D}})^2equiv 1pmod{P} ]

则此时(A^{2^{S-1}·D})又变成了在(A^2equiv1 pmod{P})的新的(A)值。
不妨设(A)值数组为({A_0,A_2,A_3,A_4....A_{S-1}})，其中(A_i=A^{2^i·D})。
（注：不难发现对于任意满足(0<i<S)的(i)，都有(A_i=2·A_{i-1})）
这样的话对于每个初始值(A)，我们可以知道所构成的(A)数组值。
此时我们可以考虑从(A_0)走到(A_{S-1})（即从(A^D)到(A^{2^{S-1}·D})）。

①考虑初始情况的(A_0)（即为(A^D)），
它如果满足(A_0equiv1 pmod{P})或(A_0equiv-1 pmod{P})。
则对于任意(A_i)（(0<i< S)），都有(A_iequiv1 pmod{P})，即通过本次测试。
（注：通过测试只能说明它可能为素数，而没通过则说明它一定不是素数。）

②考虑过程中的情况：设有一满足(0<t<S)的(t)。

1. 若有(A_tequiv-1pmod{P})，则说明
   对于任意满足(t<i<S)的(i)，都有(A_iequiv1pmod{P})，即通过此次测试。
2. 若有(A_tequiv1pmod{P})，则说明({A_{t-1}}^2equiv1pmod{P})，
   而(A_{t-1})又不满足(A_{t-1}equiv-1pmod{P})或(A_{t-1}equiv-1pmod{P})，
   因为如果(A_{t-1})满足以上条件，则在上一次或(A_0)时就已经判了。
   故存在一个(A)，使得(A^2equiv1pmod{P})且(A otin{1,-1})（定理1的延伸）。
   故未通过此次测试。

• 注：若(P(P>2))为素数： 则(A^2equiv1pmod{P})只有(A)满足
  (Aequiv1pmod{P})或(Aequiv-1pmod{P})时满足（由定理1得）。
  但(A^2equiv-1pmod{P})对于(A)的取值却没有特殊性。

③考虑结尾情况：
若对于这一组(A)都没有一个(A)使得(A^2equiv1pmod{P})，
则在判(A_{S-1})时，实际上判的是({({A^{2^{S-1}·D}})}^2 otequiv1pmod{P})
即(A^{P-1} otequiv1pmod{P})，故(P)不满足费马小定理，即(P)未通过测试。

为保证算法的正确性与高效性，
我们考虑随机化许多个整数(A)来多次判断使得(P)为素数的概率尽量大。
那么我们要随机多少次才能使得(P)素数的概率能让我们接受呢？
根据大佬们的研究表明，Prime是一个RP问题（通过大量数据来算概率），
而对于每次测试的期望失败概率是比(frac{1}{4})要小的。（通过试验得到的结论）
即我们做K次成功的测试然后判错(P)的概率即为(4^{-K})。
这个失败率在代码中也是可以接受的。
（注：经过试验，当A值取({2,3,5,7,11,13,17,19,23,29})时，可以过long long级别的询问）

代码

int MillerTest(LL N){//素数测试
	if(N<=1||N==4)return 0;//特判N比较小情况
	if(N<=3)return 1;
	LL D=N-1;
	while(D%2==0)D/=2;
	for(int i=1;i<=PK;i++){
		if(Prime[i]>=N)return 1;
		int Ok=0;
		LL A=Prime[i];
		LL x=quick_Pow(A,D,N);
		if(x==1||x==N-1)
			continue;
		while(D!=N-1){
			x=quick_Mul(x,x,N);//龟速乘防爆long long
			D*=2;
			if(x==1)return 0;
			if(x==N-1){Ok=1;break;}
		}
		if(!Ok)return 0;
	}
	return 1;
}

质因数分解

引入（同上 ）

质因数分解是OI中一个十分重要的事，在数学毒瘤题中有着举足轻重的地位。
常见的质因数分解如下：

void divide(int N){
	for(int i=2;i*i<=N;i++){
		int cnt=0;
		while(N%i==0)cnt++,N/=i;
		if(cnt)printf("%d %d
",i,cnt);
	}
	if(N!=1)printf("%d 1
",N);
}

以上是一个(O(sqrt{N}))的算法，虽然不优，但在绝大多数情况下是可以的。
但是，假如(N)的范围达到了(1e18)，以上算法很明显就不行了，我们得考虑更优的算法。
引入Pollard-Rho算法。

随机化

说起Pollard-Rho算法，我们可能不太熟悉。
但提到随机化，大家应该多少都有些了解。

//随机化找N的一个因数
int i=0;do{
	i=rand()%(N-2)+2;
}while(N%i);
printf("I found %d 
",i);

我们不难发现：随机化一次找到某个因数的概率十分小。
当(N=P·Q)时（(P，Q)均为素数），成功的概率为(frac{2}{N-1})。
期望次数为(O(N))级别，（甚至不如暴力）。

但是，如果我们加上一些优化呢？

优化1

考虑每次随机化。
其实对于随机的某个数(i)，只要(gcd(i,N) ot=1)就可以找到一个因数。
（时间复杂度仿佛还是(O(N))....）

优化2(生日悖论)

我们考虑用生日悖论来优化随机化。
尽管听起来特别玄学，但是此做法还是有一定科学依据的。

生日悖论：

生日悖论，指如果一个房间里有23个或23个以上的人，那么至少有两个人的生日相同的概率要大于50%。这就意味着在一个典型的标准小学班级(30人)中，存在两人生日相同的可能性更高。对于60或者更多的人，这种概率要大于99%。

以上是生日悖论的一个例子，其函数图像大致如下。

考虑换一种思考方式，我们随机化(N)个数出来，得到(N^2)个差值。
然后我们用这些差值来执行优化1，时间复杂度将明显降低。

考虑生日相同的情况，其可以被描述为某两个人的生日值差值为0。
从以上描述差值为0的图像中，可以明显发现选取的人数(N)越大，则概率越大。
则可以发现(N)的升大对于所有(K)的概率均有提升。
（易发现(N^2)对于每个(K)的期望值为总值域(W)的根号(O(sqrt{W}))）

然后，我们使(N)值固定。
我们设(P(K))表示以(K)为差值的概率。
则感性发现：(P(1)>P(0),P(W-1)>P(W))
感性推出：实际的(P)图像其实是在前方激增，中间接近1，后方再次猛降。
故我们在判断时，得特判差值出现在前面的情况。
（注：(N)越大时，(P)图像的前方激增斜率便越大。）

Pollard-Rho

应用优化1与优化2，考虑如何找到(N^2)个差值，
（可以认为生日悖论中的某个差值(K)是我们所需要的答案）
我们考虑随机出两个数(X)，(Y)，然后用(abs(X-Y))作为差值。
然后令(f(X)=(X^2+P) mod W,X=f(X),Y=f(Y))。
这样的话，我们只用随机化两个初值，一个常数(P)，就可以了。
但是易发现这样可能重复，即走进一个无解的循环里。
所以我们让(X,Y)的初值相等，每次随机的时候令(Y=f(f(Y)))，
即出发点一样，(Y)的速度为(X)的两倍，若有无解循环，
则总有一次(X,Y)会相遇（复杂度至多为经过的点数量），此时我们就可以直接退出了。
（注：(f(X)=(X^2+P) mod W)是一个先人传下的神奇函数）
（其神奇之处在于玄学般地大幅提升了代码效率）

其实随机这些差值的本质还是随机化，依旧是(N^2)个差值，但复杂度却玄学降低了。
???
故这个算法的复杂度十分看RP，
在一些数据的复杂度甚至超过了(O(sqrt{W})），
不过能解决一些long long的数据还是很强大了。

质因数分解

结合Miller-Rabin与Pollard-Rho算法。
我们对于一个整数(N)，先用Miller-Rabin判断其是否为质数。
若是，则上传N。
若不是，则用Pollard-Rho找一个因数出来，
将其分解成更小的两个数继续递推。

代码

//该代码只是找到了N的一个因数,可能并不是素数
LL PollardRho(LL N){//找因数
	if(N==1)return -1;
	for(int i=1;i<=PK;i++)
		if(N%Prime[i]==0)
			return Prime[i];//特判出奇迹
	LL x=(rand()%(N-2))+2,y=x;
	LL P=(rand()%(N-1))+1,D=1;
	while(D==1){
		x=(quick_Mul(x,x,N)+P)%N;
		y=(quick_Mul(y,y,N)+P)%N,y=(quick_Mul(y,y,N)+P)%N;
		D=gcd(abs(x-y),N);//x=y时,D值为N
	}
	return D;
}

例题及代码

板题：Prime Test
大意：找到(N)的最小的质因数。

#include<ctime>
#include<cstdio>
#include<algorithm>
using namespace std;
#define LL long long
const int MAXK=15;
const long long INF=1e18;
int PK=10,Prime[MAXK]={0,2,3,5,7,11,13,17,19,23,29};
//此行数对于long long范围内的素性测试已够用了
long long N;
LL Add(LL x,LL y,LL p){
	if(x+y<p)return x+y;
	return x+y-p;
}
LL gcd(LL x,LL y){//最大公因数
	if(x==0)return y;
	return gcd(y%x,x);
}
LL quick_Mul(LL x,LL y,LL p){//龟速乘
	x%=p;y%=p;
	LL ret=0;
	while(y){
		if(y&1)ret=Add(ret,x,p);
		x=Add(x,x,p);
		y>>=1;
	}
	return ret;
}
LL quick_Pow(LL x,LL y,LL p){//快速幂
	LL ret=1;
	while(y){
		if(y&1)ret=quick_Mul(ret,x,p);
		x=quick_Mul(x,x,p);
		y>>=1;
	}
	return ret;
}
int MillerTest(LL N){//素数测试
	if(N<=1||N==4)return 0;
	if(N<=3)return 1;
	LL D=N-1;
	while(D%2==0)D/=2;
	for(int i=1;i<=PK;i++){
		if(Prime[i]>=N)return 1;
		int Ok=0;
		LL A=Prime[i];
		LL x=quick_Pow(A,D,N);
		if(x==1||x==N-1)
			continue;
		while(D!=N-1){
			x=quick_Mul(x,x,N);
			D*=2;
			if(x==1)return 0;
			if(x==N-1){Ok=1;break;}
		}
		if(!Ok)return 0;
	}
	return 1;
}
LL PollardRho(LL N){//找因数
	if(N==1)return -1;
	for(int i=1;i<=PK;i++)
		if(N%Prime[i]==0)
			return Prime[i];//特判出奇迹
	LL x=(rand()%(N-2))+2,y=x;
	LL P=(rand()%(N-1))+1,D=1;
	while(D==1){
		x=(quick_Mul(x,x,N)+P)%N;
		y=(quick_Mul(y,y,N)+P)%N,y=(quick_Mul(y,y,N)+P)%N;
		D=gcd(abs(x-y),N);
	}
	return D;
}
LL Min_Div(LL N){//最小质因数
	if(N==1)return INF;
	if(MillerTest(N))return N;
	LL X=PollardRho(N);
	return min(Min_Div(X),Min_Div(N/X));
}
int main(){
	//srand(time(NULL));
	while(~scanf("%lld",&N)){
		printf("%lld
",Min_Div(N));
	}
}