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"
样例输入

样例输出

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