Miller_Rabin算法详解

提示：Miller-Rabin质数测试

小Hi：这种质数算法是基于费马小定理的一个扩展，首先我们要知道什么是费马小定理：

费马小定理：对于质数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质数测试算法。

与Fermat测试相比，增加了一个二次探测定理：

如果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测试。

但一次MR测试仍然有一定的错误率。为了使我们的结果尽可能的正确，我们需要进行多次MR测试，这样可以把错误率降低。

写成伪代码为：

Miller-Rabin(n):If (n <= 2) ThenIf (n == 2) ThenReturn TrueEnd IfReturn FalseEnd IfIf (n mod 2 == 0) Then// n为非2的偶数，直接返回合数Return FalseEnd If// 我们先找到的最小的a^u，再逐步扩大到a^(n-1)u = n - 1; // u表示指数while (u % 2 == 0) u = u / 2End While // 提取因子2For i = 1 .. S // S为设定的测试次数a = rand_Number(2, n - 1) // 随机获取一个2~n-1的数ax = a^u % nWhile (u < n) // 依次次检查每一个相邻的 a^u, a^2u, a^4u, ... a^(2^k*u)是否满足二次探测定理y = x^2 % n If (y == 1 and x != 1 and x != n - 1)// 二次探测定理// 若y = x^2 ≡ 1(mod n)// 但是 x != 1 且 x != n-1Return FalseEnd Ifx = yu = u * 2 End WhileIf (x != 1) Then// Fermat测试Return FalseEnd IfEnd ForReturn True

值得一提的是，Miller-Rabin每次测试失误的概率是1/4；进行S次后，失误的概率是4^(-S)。

题目链接：http://hihocoder.com/problemset/problem/1287

代码：

#include<cstdio>
#include<cstdlib>
#include<string>
#include<iostream>
#include<cmath>
#include<algorithm>
using namespace std;
typedef long long ll;

ll myrand()
{
    ll a = rand();
    a=a*rand();
    return a;
}

ll product_mod(ll a,ll b,ll c)    //(a*b)%c
{
    ll r=0,d=a;
    while(b>0)
    {
        if(b&1)    r=(r+d)%c;
        d=(d<<1)%c;
        b>>=1;        
    }
    return r%c;    
}

ll PowerMod(ll a,ll b,ll c)    // (a^b)%c
{
    ll r=1,d=a;
    while(b>0)
    {
        if(b&1) r=product_mod(r,d,c);
        d=product_mod(d,d,c);
        b>>=1;
    }
    return r%c;
}


bool Miller_Rabin(ll n,int times)
{
    ll k=0,i,u,x,y,a;
    if(n<2)
        return false;
    if(n==2)
        return true;
    if(!(n&1))
        return false;
    
    for(u=n-1; !(u&1); (u>>=1,++k));
    
    for(i=0;i<times;i++)
    {
        a = myrand() % (n-2)+2;
        x = PowerMod(a,u,n);
        while(u < n)
        {
            y=PowerMod(x,2,n);
            if(y == 1 && x != 1 && x != n-1)
            {
                return false;
            }
            else
            {
                x=y;
                u=u*2;    
            }            
        }
        if(x!=1)
            return false;        
    }    
    return true;
}

int main()
{    
    int t;
    ll num;
    scanf("%d",&t);    
    while(t--)
    {        
        scanf("%lld",&num);
        if(Miller_Rabin(num,20))
        {
            cout<<"Yes"<<endl;
        }
        else
        {
            cout<<"No"<<endl;
        }
    }    
    return 0;
}