Miller-Rabin质数测试

本文主要讨论使用Miller-Rabin算法编写素数的判定算法,题目来源于hihocoder

题目

题目要求

时间限制:10000ms
单点时限:1000ms
内存限制:256MB
描述
使用Miller-Rabin算法进行质数素数测试,要求输入一个数字,对其是否是素数进行判定,并打印出相对应的结果。

提示:Miller-Rabin质数测试

输入
第1行:1个正整数t,表示数字的个数,10≤t≤50
第2..t+1行:每行1个正整数,第i+1行表示正整数a[i]2≤a[i]≤10^18
输出
第1..t行:每行1个字符串,若a[i]为质数,第i行输出”Yes”,否则输出”No”
样例输入

1
2
3
4
3  
3
7
9

样例输出

1
2
3
Yes
Yes
No

题目分析

Miller-Rabin算法是一种基于费马小定理的扩展算法,首先我们需要知道什么是费马小定理,然后还要知道整个Miller-Rabin算法是如何扩展出来的。

费马小定理

费马小定理:对于质数p和任意整数a,有a^p ≡ a(mod p)(同余)。反之,若满足a^p ≡ a(mod p)p也有很大概率为质数。
将两边同时约去一个a,则有a^(p-1) ≡ 1(mod p)

也即是说:假设我们要测试n是否为质数。我们可以随机选取一个数a,然后计算a^(n-1) mod n,如果结果不为1,我们可以100%断定n不是质数。
否则我们再随机选取一个新的数a进行测试。如此反复多次,如果每次结果都是1,我们就假定n是质数。
该测试被称为Fermat测试。需要注意的是:Fermat测试不一定是准确的,有可能出现把合数误判为质数的情况。
Miller和Rabin在Fermat测试上,建立了Miller-Rabin质数测试算法。

二次探测定理

如果p是奇素数,则 x^2 ≡ 1(mod p)的解为 x ≡ 1x ≡ p - 1(mod p)
如果a^(n-1) ≡ 1 (mod n)成立,Miller-Rabin算法不是立即找另一个a进行测试,而是看n-1是不是偶数。如果n-1是偶数,另u=(n-1)/2,并检查是否满足二次探测定理即a^u ≡ 1a^u ≡ n - 1(mod n)

举个Matrix67 Blog上的例子,假设n=341,我们选取的a=2。则第一次测试时,2^340 mod 341=1。由于340是偶数,因此我们检查2^170,得到2^170 mod 341=1,满足二次探测定理。同时由于170还是偶数,因此我们进一步检查2^85 mod 341=32。此时不满足二次探测定理,因此可以判定341不为质数。

将这两条定理合起来,也就是最常见的Miller-Rabin测试。

加强版测试验证定理

尽可能提取因子2,把n-1表示成d*2^r,如果n是一个素数,那么或者a^d mod n==1,或者存在某个i使得a^(d*2^i) mod n=n-1 (0<=i<r)则我们认为n为素数。(注意i可以等于0,这就把a^d mod n=n-1的情况统一到后面去了)

这里需要注意的是,我们将该定理作为判定条件,仍然是一个不确定的概率判定条件。Miller-Rabin素性测试同样是不确定算法,我们把可以通过以a为底的Miller-Rabin测试的合数称作以a为底的强伪素数(strong pseudoprime)。第一个以2为底的强伪素数为2047。第一个以2和3为底的强伪素数则大到1 373 653。

所以我们在实际使用过程中,使用rand()函数生成随机数,或者进行多次检测判定,还是能够得到比较高的判定成功率,Miller-Rabin算法对于素数的研究判定有着巨大的辅助作用。

代码

整体代码

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
#include <iostream>
#include <cstdlib>
using namespace std;
typedef long long llong;
//求取(x * y) % n
llong mod(llong x, llong y,llong n)
{
llong res = 0;
llong temp = x % n;
while(y)
{
if(y & 0x1)
if((res += temp) > n)
res -= n;
if((temp <<= 1) > n)
temp -= n;
y >>= 1;
}
return res;
}

//求取(x ^ y) % n
llong get_mod(llong x, llong y, llong n)
{
llong res = 1;
llong temp = x;
while(y)
{
if(y & 0x1)
res = mod(res, temp, n);
temp = mod(temp, temp, n);
y >>= 1;
}
return res;
}

//编写bool函数,判定是否为素数
bool is_prime(llong n, int t)
{
if(n < 2)
return false;
if(n == 2)
return true;
if(!(n & 0x1))
return false;
llong k = 0, m, a, i;
for(m = n -1; !(m & 0x1); m >>= 1, ++k);
while(t--)
{
a = get_mod(rand() % (n - 2) + 2, m, n);
if(a != 1)
{
for(i = 0; i < k && a != n-1; ++i)
{
cout << a << endl;
a = mod(a, a, n);
}
//根据二次探测定理,只要不满足(a == 1) || (a == n - 1),就会一直遍历下去,直到最后返回false
if(i >= k)
return false;
}
}
return true;
}

//主函数
int main()
{
int times;
llong num;
cin >> times;
while(times--)
{
cin >> num;
if(is_prime(num, 1))
cout << "Yes" << endl;
else
cout << "No" << endl;
}
return 0;
}

代码分解

mod()函数

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
//求取(x * y) % n
llong mod(llong x, llong y,llong n)
{
llong res = 0;
llong temp = x % n;
while(y)
{
if(y & 0x1)
if((res += temp) > n)
res -= n;
if((temp <<= 1) > n)
temp -= n;
y >>= 1;
}
return res;
}

这个函数使用移位运算,通过将y转换成二进制形式,十分高效地求取了两个数字乘积的余数。

get_mod()函数

1
2
3
4
5
6
7
8
9
10
11
12
13
14
//求取(x ^ y) % n
llong get_mod(llong x, llong y, llong n)
{
llong res = 1;
llong temp = x;
while(y)
{
if(y & 0x1)
res = mod(res, temp, n);
temp = mod(temp, temp, n);
y >>= 1;
}
return res;
}

这个函数是经典的高次幂函数求余算法,即蒙哥马利算法,在上一篇博文中也有过介绍,博文链接
其核心思想就是将幂指数转换成二进制,通过移位运算快速地求取余数,避免了数据溢出,而且效率非常高。

is_prime()函数

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
//编写bool函数,判定是否为素数
bool is_prime(llong n, int t)
{
if(n < 2)
return false;
if(n == 2)
return true;
if(!(n & 0x1))
return false;
llong k = 0, m, a, i;
for(m = n -1; !(m & 0x1); m >>= 1, ++k);
while(t--)
{
a = get_mod(rand() % (n - 2) + 2, m, n);
if(a != 1)
{
for(i = 0; i < k && a != n-1; ++i)
{
cout << a << endl;
a = mod(a, a, n);
}
//根据二次探测定理,只要不满足(a == 1) || (a == n - 1),就会一直遍历下去,直到最后返回false
if(i >= k)
return false;
}
}
return true;
}

即数字是否是素数的判定函数,依照我们在上文提出的加强定理,包含如下要点:

  • 对所需判定的奇数n进行n-1提取因子2,把n-1表示成d*2^r的形式;
  • 取随机数a=rand(),如果a^d mod n == 1则判定为素数;
  • 如果a^d mod n != 1,则通过循环查找是否有i满足a^(d*2^i) mod n = n-1,若有,则判定为素数;
  • 如果上述条件都不成立,则遍历结果得到i == k,此时返回false

  • Tips:这里需要注意的是,Miller-Rabin算法是一个不确定算法,仍有一定的错误概率,正如上文所述的,第一个以2为底的强伪素数为2047。第一个以2和3为底的强伪素数则大到1 373 653。在一定的使用范围内仍然可以得到高效、准确的结果!
-------------本文结束感谢您的阅读-------------
0%