隐藏
Miller-Rabin与Pollard-Rho学习笔记 | Bill Yang's Blog

路终会有尽头，但视野总能看到更远的地方。

0%

Miller-Rabin与Pollard-Rho学习笔记

发表于 更新于 分类于 Valine：

Miller-Rabin

Miller-Rabin算法基于费马小定理。

费马小定理

假如$n$是素数，且$\gcd(a,n)=1$，那么$a^{n-1}\equiv1\pmod n$。

费马测试

因此，我们可以根据费马小定理得出一种检验素数的思路：
随机在$[2,n-1]$中选择一个数$a$，检验$a^{n-1}\bmod n$是否为$1$。
遗憾的是，费马小定理的逆定理不成立，也就是说，满足$a^{n-1}\equiv1\pmod n$，$n$不一定是素数。

这一类数被称为伪素数，如341能整除$2^{340}-1$，但$341=11\times31$并不是一个素数。

虽然费马小定理不成立，但伪素数并不是对所有底$a$都使得费马小定理不成立的。如上述$340$当底数$a=3$时就可以用费马小定理检验了。
因此，我们可以通过多次随机选择底数进行判断。 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
LL Quick_Pow(LL a,LL b,LL p) {
	LL sum=1;
	for(; b; b>>=1,a=a*a%p)if(b&1)sum=sum*a%p;
	return sum;
}

const int TIMES=20; 

mt19937 g(rand());

bool Miller_Rabin(LL n) {
	srand(time(NULL));
	for(int i=1; i<=TIMES; i++) {
		LL a=g()%(n-1)+1;
		if(Quick_Pow(a,n-1,n)!=1)return false;
	}
	return true;
}

Miller-Rabin素性测试

虽然上述算法看似没有太大问题，但事实上有一类数可以让其错误的概率大大提升。

卡迈克尔数（Carmichael Number）

对于合数$n$，如果对于所有正整数$b$，$b$和$n$互素，都有同余式$b^{n-1}\equiv1\pmod n$成立，则合数$n$为卡迈克尔数。

如$561=3\times11\times17$。

但我们欣喜地发现，并不是所有数都与$n$互素，如果随机选出的数不与$n$互素，即可判断出卡迈克尔数（前提是不要写成$a^p\equiv a\pmod p$）。
但，出错的概率依然很大。
有没有更精确地方法可以消除卡迈克尔数的影响？

卡迈克尔数的性质

  1. 卡迈克尔数必定是至少三个素数的乘积。
  2. 卡迈克尔数不含有素数的平方因子。

二次探测定理

对于素数$p$，且$0\lt x\lt p,e\ge1$,有：

只有两个解$x=\pm1$。

Miller-Rabin素性测试

根据卡迈克尔数的性质，可知其一定不是$p^e$。
不妨将费马小定理和二次探测定理结合起来使用：
将$n-1$分解为$n-1=u\times2^t$，不断地对$u$进行平方操作，若发现非平凡平方根时即可判断出其不是素数。

现在就能通过单次测试判断出卡迈克尔数非素数了。 

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
LL Quick_Pow(LL a,LL b,LL p) {
	LL sum=1;
	for(; b; b>>=1,a=a*a%p)if(b&1)sum=sum*a%p;
	return sum;
}

const int TIMES=10; 

mt19937 g(rand());

bool Witness(LL a,LL n) {
	LL u=n-1,t=0;
	while(u%2==0)t++,u>>=1;
	LL x=Quick_Pow(a,u,n);
	if(x==1)return false;
	for(int i=1; i<=t; i++,x=x*x%n)
		if(x!=n-1&&x!=1&&x*x%n==1)return true;
	return x!=1;
}

bool Miller_Rabin(LL n) {
	if(n==2)return true;
	if(n<2||!(n&1))return false;
	srand(time(NULL));
	for(int i=1; i<=TIMES; i++)
		if(Witness(g()%(n-1)+1,n))return false;
	return true;
}

Pollard-Rho

Pollard-Rho是用来进行大数分解的方法。
对于大数分解，我们可以进行试除，时间复杂度是$O(\sqrt n)$的，复杂度很高。

我们有一个很妙的思路，就是找到一个因子（不一定是素因子），再根据这个因子分解原数。

但是，如何找到当前数的一个因子呢？

生日悖论

$n$个人中有两个生日相同的人的概率有多大？
第$i$个人的生日有$365-i+1$种选择，故所有人的生日都不相同的概率是：

那么$n$个人中至少有两个人生日相同的概率是：

当$n=23$时，概率为$0.507$，当$n\ge60$，概率$\gt0.99$。

这看上去是荒谬的（与实际直觉相悖），因此被称为生日悖论。
生日悖论甚至被运用到密码攻击中（生日攻击）。

生日悖论的推广

对于一个数$v$，随机选一个数命中的概率为$\frac1n$，而随机选$m$个数，使它们任意两个的差命中$v$的概率约在$m=\sqrt n$时达到$\frac12$。

我们不考虑绝对值，因为绝对值只会导致碰撞的概率增大。

下图是使用Mathematica绘制的一个散点图。

因此，我们发现若选多个数，使它们的差命中$v$的概率增长速度很快。

生日悖论即可看做使差为$0$。

Pollard-Rho大数分解算法

我们尝试找到原数$n$的一个因子，方法是利用生日悖论，随机生成两个数$x,y$，若它们的差$gcd(x-y,n)\neq1$，就找到了一个因子。
根据生日悖论，我们需要将$x,y$都存下来。其实并不需要，因为$x,y$均为随机数，我们只需要找到一个随机算法生成它们即可。
这里是发明人构造的随机算法：
随机选取$x$，令其下一个数为$x^2+c$，$c$是一个随机数。
我们令$y$初始$=x$，其下一个数为$(y^2+c)^2+c$，$y$可以理解为以二倍速度行进的$x$。

$x,y$一定会循环，形成一个希腊字母$\rho$形状。

如果$x=y$，说明$y$已经遍历了环，寻找约数失败，退出重新选择$c$。
此处结合了Floyd判环法。

这样我们期望在$O(n^{\frac14})$的时间内求出$n$的一个因子。
接下来再递归因数即可得出素因子。 

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
const int TIMES=10;

LL Quick_Mul(LL a,LL b,LL p) {
	LL tmp=(a*b-(LL)((long double)a/p*b+1e-8)*p);
	if(tmp<0)return tmp+p;
	else return tmp;
}

LL Quick_Pow(LL a,LL b,LL p) {
	LL sum=1;
	for(a%=p; b; b>>=1,a=Quick_Mul(a,a,p))if(b&1)sum=Quick_Mul(sum,a,p);
	return sum;
}

mt19937 g(rand());

bool Witness(LL a,LL n) {
	LL u=n-1,t=0;
	while(u%2==0)t++,u>>=1;
	LL x=Quick_Pow(a,u,n);
	if(x==1)return false;
	for(int i=1; i<=t; i++,x=Quick_Mul(x,x,n))
		if(x!=n-1&&x!=1&&Quick_Mul(x,x,n)==1)return true;
	return x!=1;
}

bool Miller_Rabin(LL n) {
	if(n==2)return true;
	if(n<2||!(n&1))return false;
	for(int i=1; i<=TIMES; i++)
		if(Witness(g()%(n-1)+1,n))return false;
	return true;
}

LL Pollar_Rho(LL n) {
	if(!(n&1))return 2;
	while(true) {
		LL a=g()%(n-1)+1,b=a,c=g()%(n-1)+1;
		if(c==2)c=3;
		while(true) {
			a=(Quick_Mul(a,a,n)+c)%n;
			b=(Quick_Mul(b,b,n)+c)%n;
			b=(Quick_Mul(b,b,n)+c)%n;
			if(a==b)break; //环
			LL d=__gcd(abs(b-a),n);
			if(d>1)return d;
		}
	}
}

LL Find_Fac(LL n) { //找出最小因子
	if(Miller_Rabin(n))return n;
	LL p=Pollar_Rho(n);
	return min(Find_Fac(p),Find_Fac(n/p));
}

Pollard-Rho算法也有可能永远跑不出一个解，这就是玄学了（当然概率极低）。

参考资料

相关文章
姥爷们赏瓶冰阔落吧~