poj 1811 Prime Test Pollard rho算法 分解素因数 快速幂

缘起

【1】中，我们写了一个算法求出一个正整数m的欧拉函数的同时，筛出了[1,m-1]中所有与m互质的正整数, 以及输出了m分解素因数的结果（真是全能的算法诶~） 但是那是对于m较小（1e6数量级）的情况效率尚可. 但是对于超大整数，使用该算法效率就堪忧了. 于是我学习 Pollard rho 分解超大整数（long long级别）的算法.

分析

大整数分解现在任然是个世界级的难题，但却依旧是个非常重要的研究方向，这是因为大整数在公共密钥的研究上有着重要的作用。

而对于大整数的分解有很多种算法，性能上各有优异，比如大整数分解的Fermat方法，Pollard rho方法，试除法（就是我们常用的sqrt(n)那种），以及椭圆曲线法，连分数法，二次筛选法，数域分析法等等。这里，我主要讲解其中两种算法（fermat和Pollard rho）的原理和操作。

1
2
3
4
5
6
7
8
9
10
Fermat算法
首先，对于任意的一个偶数，我们都可以提取出一个2的质因子，如果结果仍为偶数，则可继续该操作，直至将其化为一个奇
数和2的多少次幂的乘积，那么我们可以假定这个奇数可以被表示成2*N+1，如果这个数是合数，那么一定可以写成N=c*d的
形式，不难发现，式中的c和d都是奇数，不妨设c>d，我们可以令a=(c+d)/2,b=(c-d)/2，那么的可以得到N=a*a-b*b，而
这正是Fermat整数分解的基础；由不等式的关系，我们又可以得到a>=sqrt(c*d)=sqrt(N)，那么，我们就可以枚举大于N
的完全平方数a*a（它增长的很快），计算a*a-N的值，判断计算的结果是否为一个完全平方数，如果是，那么a，b都是N的因子，我们就可以
将算法分治递归的进行下去，直至求出N的所有质因子。

容易看出，Fermat分解大数的效率其实并不高，但是比起试除法要好了很多；而且每次的计算都是计算出N的一个因子，更
加降低了其效率。这就让我们想着去尝试新的算法，那就是Pollard rho算法。

主要参考的是网上一篇有趣的文章《A Quick Tutorial on Pollard’s Rho Algorithm》（【2】）.

Pollard’s Rho算法最早出自 John M. Pollard 在 1975 年所写的论文

我们首先规定一下我们考虑的问题.

1
2
我们假设N是一个能被分解为 p∗q 的数
我们的目标是找到 N 的其中一个因子 p 或 q (另外一个因子可以通过除 N 来得到 )

我们首先来看看最为笨拙的试除法

1
2
3
1. 如果N是偶数，返回2. 程序结束.
2. 遍历[3,...,(int)(sqrt(N)+0.5)], 找到第一个整除N的数, 它一定就是N的最小素因子, 返回它，程序结束
3. 打印"N是素数", 程序结束.

显然试除法的复杂度是O(sqrt(N)). N超大的时候，试除法效率不高.

我们给出一个看似更加疯狂的算法

1
每次随机生成一个正整数p, 然后尝试让N去除, 如果能整除的话, 返回p, 程序结束.

你可能会说: 这特么不就是买彩票么? 这算法最坏复杂度是无限啊~ 别急，pollard rho 算法就是从lucky算法出发化腐朽为神奇的.

首先，我们定量来看看lucky算法有多糟糕.

1
2
在[2,...,N-1]这N-2个数中，假设N只有2个素因子，则lucky算法每次成功的概率是 2/(N-2)
如果N是1e10量级，则概率低于中奖.

所以我们要想办法提高命中N的质因子的概率. 学过概率论的朋友都知道所谓的生日”悖论”. 也就是一个班级中（>=60人）则几乎一定有两个人同一天过生日. 之所以叫”悖论”并不是真的是悖论，而是结论出乎我们的常识. 我们就打算使用这种birthday trick来提高命中素因子的概率.

前面之所以概率那么低是因为我们执意要直接pick N的素因子. 那当然概率低啦~ 但是如果我们改换思路——一次pick多个数，例如k个数, x1,…,xk, 判断是否存在 gcd(xi-xj, N)>1，如果存在，我们就找到了N的一个因子.

数学家证明了——约莫k=N^(1/4) 就可以保证极高的概率存在这对(xi,xj). 但是问题来了，对于超大的N（2^63 long long），即便是N^(1/4)，对内存也是一种极大的挑战. 所以Pollard rho 算法只是将2个数放进内存，其实更确切的说Pollard rho 算法并不随机生成 k 个数一股脑的放进内存并两两进行gcd(xi-xj,N)操作，而是一个一个地（随机）生成并检查连续的两个数（即这次用x生成了y, 然后求gcd(x-y,N)操作，然后用y生成z, 然后进行gcd(y-z,N)操作， 那么问题来了, 怎么随机生成呢? 可以使用一个函数 $f(x)=x^2+a \ (mod \ N)$ 来根据x生成y）。 反复执行这个步骤并希望能够得到某次求gcd的结果能>1（因为gcd算法是欧几里得，效率很快）

但是这里有个问题，万一随机生成（例如，使用上面的$f(x)=x^2+a \ (mod \ N)$随机生成随机数）的随机数陷入了循环怎么办? 例如我们构造出来的随机数为

2, 10, 16, 23, 29, 13, 16, 23, 29, 13 … …
在这个例子中，我们最终将会在16, 23, 29, 13这个圈中无限循环

虽然大部分情况不会出现这种循环, 但只要我们无法严格证明不会出现环，就必须保证程序的健壮性.

1
2
ps: 为什么叫 pollard rho 算法? pollard自然不消说——该算法的发明人就是pollard. rho 就是希腊字母, 
就是一个环状. 而且，此算法最得意之处就是环状的处理上（floyd的想法）

也就是问: 如何探测环的出现呢？ 显然，一种笨办法是记录所有出现过的数，一旦出现了之前出现过的数，只要产生随机数的函数没有变化，则随机序列一定出现环. 这种办法是不切实际的（内存不够）. 正确的方法是Floyd发明的一个算法，它是一个极为巧妙和聪明的方法. 假设我们在一个周长很长的跑道上跑步，你怎么知道你已经跑完一圈了呢? 有人会说——我发现我第二次经过一个旗杆了，嗯，这种方法就是刚刚说的记录所有出现过的数来判断是否出现环. 更加聪明的办法是找一个速度恰好是我2倍的人和我同方向，同时出发，一旦我被套圈，则说明套圈的那一刹那，我一定恰好跑完了一圈. Floyd只需要用$f$和$f(f)$ （后者就是前者速度的2倍）同时发生随机数. 一旦出现生成的随机数相等，则一定出现了环. 则pollard rho算法要退出（因为已经陷入死循环了）并且重新选择随机函数f.

好了，pollard rho 算法的思想已经讲完了， 我们来看看 poj 1811 Prime Test

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
给定一个超大的正整数，你需要确定它是不是素数

【输入】
多样例，第一行是样例数T， 后面每一行都是一个N (2 <= N < 2^54)

【输出】
对每个样例，如果N是素数，输出"Prime"， 如果不素数，输出N的最小素因子

【样例输入】
2
5
10

【样例输出】
Prime
2

【限制】
Time limit 6000 ms
Case time limit 4000 ms

判定素数显然使用 Miller-Rabin测试. 但是注意一点，这里N极大，不能像【3】中那样只使用 2、7、61三个素数就可以进行素性测试了（因为【3】中N的范围仅仅到10亿）

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
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
//#include "stdafx.h"

#include <cstdio>
#include <cstdlib>
#include <ctime>
#include <cmath>
//#define LOCAL

typedef long long LL;
const int test = 12, A = 240; // 进行Miller-Rabin测试12次, A就是pollard rho算法中的随机生成函数f中的a(注意，这里不可随机，不然吃T)
LL ret,n; // n的最小素因子是ret

LL modu_mul(LL a, LL b, LL mod) // a*b 模掉mod 的余数
{
	a %= mod, b%=mod;
	LL ans = 0;
	while(b)
	{
		if (b&1)
		{
			ans = (ans+a)%mod;
		}
		if (b>1)
		{
			a = (a<<1)%mod;
		}
		b>>=1;
	}
	return ans;
}

LL modu_pow(LL a, LL b, LL mod) // a^b模掉mod的余数
{
	LL ans = 1;
	while(b)
	{
		if (b&1)
		{
			ans = modu_mul(ans, a, mod);
		}
		if (b>1)
		{
			a = modu_mul(a, a, mod);
		}
		b>>=1;
	}
	return ans;
}

bool mr(LL p)
{
	if (p<2)
	{
		return false;
	}
	if (!(p&1))
	{
		return p==2;
	}
	LL s = 0, t = p-1;
	while(!(t&1))
	{
		s++;
		t>>=1;
	} // p-1=t*2^s
	for (int i = 0; i<test; i++)
	{
		LL a = rand()%(p-1)+1; // [1,p-1]中随机一个数
		LL x = modu_pow(a,t,p); // a^t mod p =x
		for (int j = 0; j<s; j++)
		{
			LL y = modu_mul(x,x,p);
			if (y==1 && x!=1 && x!=p-1) // Miller-Rabin测试
			{
				return false;
			}
			x = y;
		}
		if (x!=1) // 费马素性测试
		{
			return false;
		}
	} // 通过了 test次Miller-Rabin测试
	return true;
}

LL gcd(LL a, LL b)
{
	return b?gcd(b,a%b):a;
}

LL pollard_rho(LL n, LL a) // pollard rho算法 其中a是文章中的那个随机数发生函数f中的a, 这里对n采用pollard rho算法, 返回n的一个非平凡因子, 但是如果算法失败（即随机发生函数f出现环），则返回n
{
	LL i = 1, k=2, x = rand()%n, y =x;
	while (1)
	{
		i++;
		x = (modu_mul(x,x,n)+a)%n; // f(x)=x^2+a mod n
		LL d = gcd(y-x+n, n); // pollard rho 的增大概率的想法,这里没用abs(x-y),因为效率低下
		if (1<d&&d<n)
		{
			return d; // 返回非平凡因子
		}
		if (x==y) // 出现环了（如果有环的话, 则一定会发生x==y的，因为x的速度是y的2倍）
		{
			return n;
		}
		if (i==k)
		{
			y=x;
			k<<=1; // 如果我们把进行一次f运算视作走了一格的话，则x走的速度是y的两倍
		}
	}
}

void sz(LL n, LL a) // 递归寻找最小素因子
{
	if (mr(n)) // 如果n是素数（其实这里就可以收集n的素因子了）
	{
		if (n<ret)
		{
			ret = n;
		}
		return;
	}
	LL d = n; // n是合数， 则就一定能从n中通过pollard rho算法榨出一个非平凡因子来
	while(d==n)
	{
		d = pollard_rho(n, a--); // 如果pollard rho算法失败, 则算法返回n, 则就要重新选择pollard rho中的随机生成函数f中的a（这里是直接减1）
	}
	sz(d, a);
	sz(n/d, a); // 递归去改变ans
}

int main() 
{
#ifdef LOCAL
	freopen("d:\\data.in", "r", stdin);
	freopen("d:\\my.out", "w", stdout);
#endif
	srand((int)(time(0)));
	int t;
	scanf("%d", &t);
	while(t--)
	{
		scanf("%lld", &n);
		if (mr(n))
		{
			puts("Prime");
		}
		else
		{
			if (!(n&1))
			{
				puts("2");
			}
			else if ((n&1)&!(n%3))
			{
				puts("3");
			}
			else
			{
				ret = n; // 初始化最小素因子为n
				sz(n,A);
				printf("%lld\n",ret);
			}
		}
	}
	return 0;
}

ac情况

Status Accepted
Time 969ms
Memory 88kB
Length 2479
Lang C++
Submitted 2019-08-17 22:32:39
Shared
RemoteRunId 20759716

参考

【1】https://yfsyfs.github.io/2019/08/16/poj-2773-Happy-2006-%E7%BA%BF%E6%80%A7%E7%AD%9B/#more

【2】http://www.cs.colorado.edu/~srirams/classes/doku.php/pollard_rho_tutorial

【3】https://yfsyfs.github.io/2019/08/16/POJ-3641-Pseudoprime-numbers-Miller-Rabin%E7%B4%A0%E6%95%B0%E6%B5%8B%E8%AF%95/