MR素性探测
简介
MR算法全称是Miller-Rabin测试,是一个非确定的算法,用于判断一个数是否是质数.虽然是一个非确定的算法,但是只要巧妙地选取参数,在一定范围内就是一个确定性的算法.
前置条件:费马小定理 1 ≡ a^(p-1) (mod p)
Miller和Rabin两个人的工作让Fermat素性测试迈出了革命性的一步，建立了传说中的Miller-Rabin素性测试算法。 新的测试基于下面的定理(
二次探测定理)：如果p是素数，x是小于p的正整数，且x^2 mod p = 1，那么要么x=1，要么x=p-1。这是显然的，因为x^2 mod p = 1相当于p能整除x^2-1，也即p能整除(x+1)(x-1)。由于p是素数，那么只可能是x-1能被p整除(此时x=1)或x+1能被p整除(此时 x=p-1)。
这就 是Miller-Rabin素性测试的方法。不断地提取指数n-1中的因子2，把n-1表示成
d*2^r（其中d是一个奇数）。那么我们需要计算的东西就 变成了a的d*2^r次方除以n的余数。于是，a^(d * 2^(r-1))要么等于1，要么等于n-1。如果a^(d * 2^(r-1))等于1，定理继续适用于a^(d * 2^(r-2))，这样不断开方开下去，直到对于某个i满足a^(d * 2^i) mod n = n-1或者最后指数中的2用完了得到的a^d mod n=1或n-1。这样，Fermat小定理加强为如下形式：尽可能提取因子2， 把n-1表示成d*2^r，如果n是一个素数，那么或者a^d mod n=1，或者存在某个i使得a^(d*2^i) mod n=n-1 ( 0<=i<r )（注意i可以等于0，这就把a^d mod n=n-1的情况统一到后面去了）
关于正确性
如果你每次都用前7个素数(
2, 3, 5, 7, 11, 13和17)进行测试，所有不超过341 550 071 728 320(即3.4e14)的数都是正确的。如果选用2, 3, 7, 61和24251作为底数，那么10^16内唯一的强伪素数为46 856 248 255 981(即4.6e13)。
注意是当取模为1时才能继续降幂取模,取模为N-1时算法终止(但是之后的算法实际上不是这么死板地实现的).
模板
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
ll Quick_Multiply(ll a, ll b, ll c) //快速积（和快速幂差不多）
{
long long ans = 0, res = a;
while (b) {
if (b & 1)
ans = (ans + res) % c;
res = (res + res) % c;
b >>= 1;
}
return ans;
}
ll Quick_Power(ll a, ll b, ll c) //快速幂，这里就不赘述了
{
ll ans = 1ll, res = a;
while (b) {
if (b & 1)
ans = Quick_Multiply(ans, res, c);
res = Quick_Multiply(res, res, c);
b >>= 1;
}
return ans;
}
bool Miller_Rabin(ll n, int a, ll d) //判断素数
{
if (n == 2) return true;
if (n == a) return true;
if ((n & 1) == 0) return false;
while (!(d & 1)) {
d >>= 1;
}
ll t = Quick_Power(a, d, n);
while ((d != n - 1) && (t != 1) && (t != n - 1)) {
t = (t * t) % n; // 注意此处可能会产生溢出,取决于n的大小,如果n大于1e9,估计就gg
// t = Quick_Multiply(t, t, n); // 避免了溢出,但是会导致整体复杂度加一个Log1e18,大概是60左右的常数
d <<= 1;
}
return (t == n - 1 || (d & 1) == 1);
}
bool isPrime(ll n) {
if (n < 2) return false;
int MRprime[] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29};
for (int i = 0; i < 10; ++i) {
if (!Miller_Rabin(n, MRprime[i], n - 1)) {
return false;
}
}
return true;
}
POj 3641
题意
给你一个p,一个a,问p是不是基于a的伪素数.
解析
伪素数是指根据费马小定理判断是素数的,实际上却不是素数的数.
这题的解法很显然,就是先用费马小定理判断一下是不是,再真正地判断一下是不是.而真正的判断的话,很多题解上给出的是暴力判断,我们就以此当做MR素性测试的模板题吧.
费马小定理前文提到了,但是在实际过程中发现了一点问题:
如果根据1 ≡ a^(p-1) (mod p),居然会WA,但是如果根据p ≡ a^(p) (mod p),就能过.这实在是匪夷所思,之后再深入探讨.
AC代码
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
#include <iostream>
using namespace std;
#define rd(x) (rand()%(x))
typedef long long ll;
int prime[] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29};
ll Quick_Multiply(ll a, ll b, ll c) //快速积（和快速幂差不多）
{
long long ans = 0, res = a;
while (b) {
if (b & 1)
ans = (ans + res) % c;
res = (res + res) % c;
b >>= 1;
}
return ans;
}
ll Quick_Power(ll a, ll b, ll c) //快速幂，这里就不赘述了
{
ll ans = 1ll, res = a;
while (b) {
if (b & 1)
ans = Quick_Multiply(ans, res, c);
res = Quick_Multiply(res, res, c);
b >>= 1;
}
return ans;
}
bool Miller_Rabin(ll n, int a, ll d) //判断素数
{
if (n == 2) return true;
if (n == a) return true;
if ((n & 1) == 0) return false;
while (!(d & 1)) {
d >>= 1;
}
ll t = Quick_Power(a, d, n);
// 这个循环运用了一些技巧,
// 尽可能提取因子2， 把n-1表示成d*2^r
// 如果n是一个素数，
// 那么或者a^d mod n=1，
// 或者存在某个i使得a^(d*2^i) mod n=n-1 ( 0<=i<r)
// 循环的次数应该是r,也就是初始的时候d移除后面连续的所有0,最后再循环里一次一次挪回来
// 如果d又重新变成了N-1,那么说明循环终止
// 此处的t就是a^(d*2^i) mod n ( 0<=i<r) 只要他在循环中途等于n-1,则N一定是质数
// 还有一种可能是a^d mod n == 1
// 这种情况只会发生在第一次循环中,所以如果是因为t==1而跳出循环,则需要判断d是否是奇数
// 因为只有在循环的第一轮(d为奇数),t==1才能认为是n是质数的证据
while ((d != n - 1) && (t != 1) && (t != n - 1)) {
t = (t * t) % n; // 注意此处可能会产生溢出,取决于n的大小,如果n大于1e9,估计就gg
// t = Quick_Multiply(t, t, n); // 避免了溢出,但是会导致整体复杂度加一个Log1e18,大概是60左右的常数
d <<= 1;
}
return (t == n - 1 || (d & 1) == 1);
}
bool isPrime(ll n) {
if (n < 2) return false;
int MRprime[] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29};
for (int i = 0; i < 10; ++i) {
if (!Miller_Rabin(n, MRprime[i], n - 1)) {
return false;
}
}
return true;
}
int main() {
ll N, P;
while (cin >> P >> N) {
if (N == 0 && P == 0) {
break;
}
bool f1 = isPrime(P);
bool f2 = Quick_Power(N, P, P) == N; // 此处实在是匪夷所思 写成 N^(P-1) == 1 (mod P) 就会WA,但是这样写就没问题
if (f1) {
cout << "no" << endl;
} else if (f2) {
cout << "yes" << endl;
} else {
cout << "no" << endl;
}
}
return 0;
}
