Miller-Rabin & Pollard-Rho

发表于 | 更新于

Problem

质因子分解一个给定的很大的NN

Solution

Part I:质数判定 Miller-Rabin算法

费马小定理

，逆定理成立的概率也很大。 如果任意选择一个数字a0a_0，发现 ，那么PP就有很大概率是个质数。如果随机地多取几个a0a_0，那么出错的概率会大大减小。 但是仍然有一些合数XX，满足对于 。这些数字成为Carmichael数。为了把这些数字也正确地判断掉，引入了

二次探测定理

PP是一个质数，那么 的解 。 也就是说如果 且 ，那么XX不是质数。

做法

所以就随机选几个数字a[]a[]，用费马小定理结合二次探测定理来测试一下。 具体的做法，就是提取P−1P-1中的22，把P−1P-1化为 的形式 求出ada^d，一路平方，如果出现了 而 ，就不是质数。 最后会得到aP−1a^{P-1}，再看看这个是不是=1=1。

Part II:寻找因子 Pollard-Rho算法

寻找因子，试除法的复杂度是O(N)O(\sqrt N)的，处理不了64位整数的范围。 于是Pollard-Rho诞生了

生日悖论

从[1,N][1,N]里面选两个数，两个数字相等的概率是1N2\dfrac 1{N^2} 但如果从[1,N][1,N]里面选三个数，使得前两个数字的差等于第三个数，那么概率就能增加不少，因为每个数对应着两个和它差值为一个定值的数。 如果选的更多，那么概率就更大。用一个通俗的模型解释，就是在≥23\ge 23人中，有>50%> 50\%的概率出现两个人生日在同一天，然后这似乎并不符合常理，所以就是生日悖论。

Pollard-Rho

提到了生日悖论，多半是要用基于概率的算法了。。大体的做法就是遍历一个随机数序列，看相邻两项的差与给定xx是否互质，如果不互质，就找到了一个因数。 因为随机数序列存不下，那就一项一项地生成，但是随机数序列会有周期，会卡在里面出不来。 而由于空间限制，判周期只能用Floyd套圈法，那就稍微改变一下策略，让两个指针在随机数序列上分别一次走一步、走两步，直到差与xx不互质。(注：gcd(0,x)==x\gcd(0,x)==x) 在实现的时候，有一个优化。两个指针同时走，会重复地计算很多个项。那么只让一个指针走，另一个指针在它走了2k2^k步之后一步走过去。具体见代码。这个优化的正确性应该比较显然，但是还没有细想会不会产生一些退化的情况。

Code

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
#include "lucida"
const int prime[]={2,3,5,7,11,13,17,19,23,29},C_MAX=99929;
using std::swap;
using std::max;
/*
int64 mul(int64 x,int64 times,int64 modu) {
	if(x<times) swap(x,times);
	int64 res;
	for(res=0;times;(x+=x)%=modu,times>>=1)
		if(times&1) (res+=x)%=modu;
	return res;
}
*/
int64 mul(int64 a,int64 b,int64 modu)
{
    int64 temp=a*b-(int64)((long double)a/modu*b+1e-8)*modu;
    if (temp<0) temp+=modu;
    return temp;
}

int64 pow(int64 base,int64 index,int64 modu) {
	int64 res;
	for(res=1;index;base=mul(base,base,modu),index>>=1)
		if(index&1) res=mul(res,base,modu);
	return res;
}
int64 gcd(int64 a,int64 b) {
	a=a<0?-a:a;b=b<0?-b:b;
	for(int64 t;b;t=a%b,a=b,b=t);
	return a;
}
bool isPrime(int64 x) {
	if(x==2) return 1;
	if(x<2 || ~x&1) return 0;
	int64 s,d;
	for(d=x-1,s=0;~d&1;s++,d>>=1);
	for(int i=0;i<9;++i) {
		if(x==prime[i]) 
			return 1;
		int64 cur=pow(prime[i],d,x);
		if(cur==1 || cur==x-1) continue;
		int64 pre=cur;
		for(int k=1;k<=s;++k) {
			cur=mul(cur,cur,x);
			if(cur==1 && pre!=1 && pre!=x-1) 
				return 0;
			pre=cur;
		}
		if(pre!=1) 
			return 0;
	}
	return 1;
}
int64 FindDivisor(int64 x) {
	int64 c=rand()%C_MAX,d;
	for(int64 p2=rand()%x,p1=(mul(p2,p2,x)+c)%x,k=2,s=1;(d=gcd(p1-p2,x))==1;++s,p1=(mul(p1,p1,x)+c)%x)
		if(k==s) {
			p2=p1;
			k<<=1;
		}
	return d;
}
int64 Solve(int64 x) {
	//os<<x<<'\n';
	if(isPrime(x)) return x;
	int64 r;while((r=FindDivisor(x))==x);
	return max(Solve(r),Solve(x/r));
}
void Work() {
	int64 x;is>>x;
	int64 pd=Solve(x);
	if(pd==x) os<<"Prime"<<'\n';
	else os<<pd<<'\n';
}
int main() {
	freopen("input","r",stdin);
	srand(0x1f1f1f1f);
	int T;is>>T;
	while(T--)
		Work();
	return 0;
}