Miller-Rabin质数测试
本文主要讨论使用Miller-Rabin算法编写素数的判定算法，题目来源于hihocoder。
题目
题目要求
时间限制:10000ms
单点时限:1000ms
内存限制:256MB
描述
使用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"
样例输入
3
3
7
9
样例输出
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 ≡ 1或x ≡ p - 1(mod p)
如果a^(n-1) ≡ 1 (mod n)成立，Miller-Rabin算法不是立即找另一个a进行测试，而是看n-1是不是偶数。如果n-1是偶数，另u=(n-1)/2，并检查是否满足二次探测定理即a^u ≡ 1或a^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算法对于素数的研究判定有着巨大的辅助作用。
代码
整体代码
#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()函数
//求取(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()函数
//求取(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()函数
//编写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。在一定的使用范围内仍然可以得到高效、准确的结果！
Github： https://github.com/haoyuanliu
个人博客： http://haoyuanliu.github.io/
个人站点，欢迎访问，欢迎评论！