POJ 3641 Pseudoprime numbers Miller-Rabin素数测试

缘起

素性测试(即测试给定的数是否为素数)是近代密码学中的一个非常重要的课题.

而众所周知，测试一个数是不是素数，我们有O(sqrt(N))算法. 但是对于超大的n, 这个算法复杂度依旧很高. 所以我们来学习著名的Miller-Rabin素数测试算法。这是一种基于概率的算法,说白了，就是如果一个数能全部通过次数相当的Miller-Rabin测试的话，则我们以非常高的置信度判定它是质数.

分析

首先有费马小定理与二测探测定理(为什么叫二次探测定理，看后面的Miller-Rabin测试算法就知道了)如下

1
2
3
4
5
6
7
费马小定理
p是素数, a是整数, 且(a,p)=1, 则 a^(p-1)模掉p等于1.
证明: 考虑{1,...,p-1}和{a,...,a*(p-1)}, 后者中任意两个彼此不模p同余. 所以1*2*...*(p-1)和a*2a*..*(p-1)a 模p同余. 即 (a^(p-1)-1)*(p-1)!可以整除p, 所以 a^(p-1)模p余1. 费马小定理得证.

二次探测引理
p是一个素数, x^2模p等于1的话, 则要么x模p为1要么模p为p-1
证明: x^2-1=(x+1)*(x-1), 可以被p整除. 得证.

基于费马小定理的素数测试是费马测试, 费马测试的算法如下

1
2
3
4
5
6
7
1 随机取一个 a
2 如果它不满足 a^(n-1)%n ==1
3 则n一定是 合数
4 退出
5 如果它满足 a^(n-1)%n ==1
6 则它是一个素数的概率是1/2
7 回到 1

如果对这个过程重复50次，每次都没说它是合数，那这个数是素数的概率是多少那？ 它只有(1/2)^50可能不是素数！

很显然就像是投硬币那样如果你投了1000次，每次都投到正面 ，那你就会怀疑这个硬币2面都是正面！

但是素数的费马测试算法是有漏洞的,因为有一种数叫做卡米切尔数，卡米切尔数使得上面发生成为大概率事件,例如561(卡米切尔数，561=3*11*17)有318个伪证据(即底数a),使得将561判为素数

因此费马素性测试几乎就没用了. 那既然费马小定理只是素数的必要条件，所以才出的卡米切尔数这种漏洞. 那么有没有充要条件呢? 有的！Wilson爵士给出了如下定理

1
2
Wilson定理:
对于给定的正整数p，p是素数的充要条件为 (p-1)! mod p 的结果为 p-1

Wilson定理的证明是容易的（参考【2】）

对于p=2，Wilson定理的成立是显然的，我们下面假设p是奇数

充分性: 如果p是素数，则 (p-1)! = 1*(p-1) *(2*3*…*(p-2)) , 只需要论证2*3…\(p-2)（偶数个数连乘）除以p余数为1即可. 我们其实只需要证明对于任意 $a\in {2,…,p-1}$ , 存在 $b\in {2,…,p-1}$, 且 b!=a, 使得 a*b 除以p余数为1.

只需要考虑集合 {1*a, 2*a,…., (p-1)*a}, 这(p-1)个数构成的新的集合. 首先他们两两之间不会mod p同余(因为p是素数). 但是又是p-1个数，所以和{1,…,p-1}模掉p的结果是一样的. 存在b, 使得 b*a模掉p余数为1. 但是a会等于b吗? 不会, 因为根据上面的二次探测定理 , 如果a=b的话，则a只能模掉p为1或者p-1, 则a只能是1或者p-1,但是伊始我们假设$a\in {2,…,p-1}$ 的. 所以b一定不等于a，即得证.

必要性: 必要性就简单多了. 反证法，如果p不是素数的话，p=a*b（a>1, b>1）, 则(p-1)!和p之间有公约数a和b, 从而

a和b 都是(p-1)的因子，所以(p,p-1)>1 ，但是 p和p-1显然是互质的.

Wilson定理得证.

Wilson定理给出了素数的充要条件诶~ 但遗憾的是根据Wilson定理来素性测试所需的计算量太大，无法实现对较大整数的测试. 所以Wilson定理也不能用于素性测试。于是才需要引入今天的主角——Miller-Rabin 素性测试（主要援引的是【1】）.

1
2
3
4
令 p 为一个奇素数, p-1 = t*2^s, t>1是奇数, s是正整数.
则考虑任取一个小于p的正整数a. 则考虑如下的序列（称之为Miller序列）
m={a^t, a^(2t),a^(4t),...,a^(t*2^(s-1))}, 即m[i] = a^(t*2^(i)), 0<=i<s, Miller序列共计s项
如果m[i]模掉p为1，但是m[i-1]模掉p既不是1又不是p-1的话（s>i>=1）, 则根据二次探测定理, p一定不是素数.

原理讲清楚了，下面解决poj 3641

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
费马定理指出，对于任何素数p和任何整数a> 1，a^p = a（mod p）。 p的一些（但不是很多）非素数值，称为base-a 
伪素数，对于某些a具有此属性。 （有些人称为Carmichael数字，这些数字对于所有的a都是base-a伪素数。）

给定2 <p≤1000000000且1 <a <p，确定p是否是base-a 伪素数。

【输入】
多样例，以0 0结尾，每个样例占一行，分别表示p和a

【输出】
对于每个样例，如果p是base-a伪素数，则输出yes, 否则输出no

【样例输入】
3 2
10 3
341 2
341 3
1105 2
1105 3
0 0

【样例输出】
no
no
yes
no
yes
yes

【限制】
Time limit 1000 ms
Memory limit 65536 kB

题意就是

1
2
由费马小定理可得，对于素数p，a^ p  = a (mod p)，但是对于某些非素数p,也有比较小的可能满足a^ p  = a (mod 
p)，如果满足，则称p是a条件下的伪素数,现给出p，a，问p是不是a条件的伪素数。

解法就是首先使用Miller-Rabin测试p是不是素数，如果是，则肯定输出no（因为伪素数必须不是素数），如果不是的话，再用快速幂求a^p除以p的余数, 看看是不是a, 如果是的话，则p就是base-a伪素数, 输出yes，否则输出no

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
//#include "stdafx.h"

#include <stdio.h>
//#define LOCAL
typedef long long LL;

LL p,a;
const LL primes[3] = {2,7,61}; // 三轮Miller-Rabin测试足以

LL modu_mul(LL x, LL y, LL mod) // 求x*y模掉mod的余数, 注意, 不能直接乘, 因为LL * LL会爆掉, 所以转乘法为加法（沿用快速幂的思想）
{
	LL ans = 0;
	x%=mod,y%=mod; // 能少一些是一些
	while(y) // 效仿快速幂
	{
		if (y&1)
		{
			ans = (ans+x)%mod;
		}
		if (y>1)
		{
			x = (x<<1)%mod;
		}
		y>>=1;
	}
	return ans;
}

LL modu_pow(LL x, LL y, LL mod) // 求 x^y 模掉mod的余数
{
	LL ans = 1;
	while(y) // 快速幂
	{
		if (y&1)
		{
			ans = modu_mul(ans, x, mod);
		}
		if (y>1)
		{
			x = modu_mul(x,x, mod);
		}
		y>>=1;
	}
	return ans;
}

bool mr() // Miller-Rabin 素性测试模板, 返回true表示p是素数
{
	if (p<2)
	{
		return false;
	}
	if (!(p&1))
	{
		return p==2;
	}
	LL s = 0, t = p-1; // 假设p是大于1的奇数，则考虑分解p-1=t*2^s
	while (!(t&1))
	{
		t>>=1;
		s++;
	} // 最后 p-1=t*2^s（s>=0, t是大于1的奇数）
	for (LL i = 0; i<3;i++) // 进行三轮Miller-Rabin测试
	{
		if (p==primes[i])
		{
			return true;
		}
		LL x = modu_pow(primes[i], t, p); // x 是 primes[i]^t模掉p的余数
		for (LL j = 0; j<s; j++)
		{
			LL y = modu_mul(x,x,p);
			if (y==1 && x!=1 && x!=p-1)
			{
				return false;
			}
			x = y;
		}
		if (x!=1) // 费马素数测试——如果x=prime[i]^(p-1)模p的余数不为1的话,则p一定不是素数(by 费马小定理)
		{
			return false;
		}
	}
	return true; // 通过了3次Miller-Rabin测试, 以很高置信度判定p是素数
}

int main() 
{

#ifdef LOCAL
	freopen("d:\\data.in", "r", stdin);
#endif
	while(scanf("%lld%lld", &p, &a), p&&a)
	{
		if (mr()) // 如果p是素数
		{
			puts("no");
		}
		else if(modu_pow(a,p,p)==a)
		{
			puts("yes");
		}
		else puts("no");
	}
	return 0;
}

ac情况

Status Accepted
Time 16ms
Memory 80kB
Length 1538
Lang C++
Submitted 2019-08-17 11:36:01
Shared
RemoteRunId 20756077

参考

【1】https://blog.csdn.net/ECNU_LZJ/article/details/72675595

【2】http://blog.sina.com.cn/s/blog_14a210e280102x9il.html