Miller-Rabin 素数测试

· · 题解

众所周知，判断一个数是否是素数有一个 O(\sqrt{n}) 的算法：

inline bool isPrime(ll n)
{
    if(n<=1)return false;
    for(int i=2;i<=sqrt(n);i++)
        if(n%i==0)return false;
    return true;
}

那么，有没有什么更快的算法呢？

那就是 Miller-Rabin 素数测试！

Miller-Rabin 素数测试基于这个定理（ex费马小定理）：

如果 n 为素数，则 \forall a<n，设 n-1=d\cdot 2^r，d 为奇数，则以下两个命题至少有一个成立：

  1. a^d\equiv 1\pmod n
  2. \exists 0\le i<r,\text{s.t. }a^{d\cdot 2^i}\equiv -1\pmod n

请注意，它的逆命题不一定成立。

Miller-Rabin 算法是这样的：

取一些数 1\le a_1,a_2,\cdots,a_k<n，然后把它们作为定理中的 a 代入到上述定理中看看定理是否成立，如果都成立，那么 n 是素数，否则 n 不是素数。

你可能会问：这个算法得到的结果不是会有 \dfrac{1}{2^k} 的概率是错的吗？事实证明，只要你选的数强度足够大（例如 [2,3,5,7,11,13,17,37]）那么这个算法得到的结果在 long long 范围内错误的概率是 \large{0\%}。

所以，Miller-Rabin 就是如此简单。它的时间复杂度为 O(k\log n)。

Miller-Rabin 模板题：SP288

代码：（注意要用龟速乘）

#include<cstdio>
#include<algorithm>

typedef long long ll;

const ll prime[]={2,3,5,7,11,13,17,37};

ll mul(ll x,ll y,ll p)
{
    ll res=0;
    while(y)
    {
        if(y&1)res=(res+x)%p;
        x=(x+x)%p;
        y>>=1;
    }
    return res;
}

ll qpow(ll x,ll y,ll p)
{
    ll res=1;
    while(y)
    {
        if(y&1)res=mul(res,x,p);
        x=mul(x,x,p);
        y>>=1;
    }
    return res;
}

bool MillerRabin(ll n,ll a)
{
    ll d=n-1,r=0;
    while(!(d&1))
        d>>=1,r++;
    ll z=qpow(a,d,n);
    if(z==1)return true; // 情况 1
    // 情况 2
    for(int i=0;i<r;i++)
    {
        if(z==n-1)return true;
        z=mul(z,z,n);
    }
    return false;
}

bool isPrime(ll n)
{
    if(n<=1)return false;
    for(int i=0;i<8;i++)
    {
        if(n==prime[i])return true; // 特判 1
        if(n%prime[i]==0)return false; // 特判 2
        if(!MillerRabin(n,prime[i]))return false;
    }
    return true;
}

int main()
{
    int T;
    scanf("%d",&T);
    while(T--)
    {
        ll n;
        scanf("%lld",&n);
        if(isPrime(n))puts("YES");
        else puts("NO");
    }
    return 0;
}