Miller-Rabin 素性测试算法

前提

Miller-Rabin 素性测试算法需要如下两个引理：

1. 费马小定理

设 $p$ 是素数，$a$ 为整数，且 $(a,p)=1$，则 $a^{p-1}\equiv1\pmod{p}$

证明:

$\begin{aligned} &考虑\ 1,2,3,\dots,(p-1)\ 这\ p-1\ 个数字，同时乘\ a\ 得\ a,2a,3a,\dots,(p-1)a。\\ &\because a\equiv b\pmod{p},(c,p)=1\\ &\therefore ac\equiv bc\pmod{p} \\ &\therefore 1\times2\times\dots\times(p-1)\equiv a\times2a\times\dots\times(p-1)a \pmod{p}\\ &\therefore (p-1)!\equiv (p-1)!~\times a^{p-1}\pmod{p}\\ &\because ((p-1),p)=1\\ &\therefore a^{p-1}\equiv1\pmod{p} \end{aligned}$

2. 二次探测定理

设 $p$ 是一个素数，且 $0<x<p$，则方程 $x^2\equiv1\pmod{p}$ 的解为 $x=1$ 或 $p-1$。
证明：

$\begin{aligned} &易知\ x^2-1\equiv0\pmod{p}\\ &\therefore (x+1)(x-1)\equiv 0 \pmod{p}\\ &\therefore p|(x+1)(x-1)\\ &\because p\ 是素数\\ &\therefore x=1\ 或\ x=p-1 \end{aligned}$

算法介绍

首先，费马小定理只是判定 $p$ 为质数的必要条件。即费马小定理不成立，$n$ 一定是合数；费马小定理成立，$n$ 也有可能不是质数。

下面来看 Miller-Rabin 算法的分析过程：

假设 $n$ 为素数且 $n>2$，则 $n-1$ 为偶数。令 $n-1=2^q\cdot m$。

随机选取整数 $a(0<a<n)$。

由费马小定理可得 $(a^{2^q\cdot m}=a^{n-1})\equiv 1\pmod{n}$。

再由二次探测定理可得 $a^{2^{q-1}\cdot m}\equiv 1\pmod{n}$ 或 $a^{2^{q-1}\cdot m}\equiv n-1\pmod{n}$。

若 $a^{2^{q-1}\cdot m}\equiv 1\pmod{n}$ 成立，那么再由二次探测定理可得 $a^{2^{q-2}\cdot m}\equiv 1\pmod{n}$ 或 $a^{2^{q-2}\cdot m}\equiv n-1\pmod{n}$……

反复使用二次探测定理拆解，直到拆解到 $a^m\equiv 1\pmod{n}$ 或 $a^m\equiv n-1\pmod{n}$。

总结一下：$a^m\equiv 1\pmod{n}$ 或存在 $0\le r<q$ 使得 $a^{2^r\cdot m}\equiv n-1\pmod{n}$，则称 $n$ 通过 Miller 测试。

但是还存在一种伪素数，它们可以通过 Miller 测试但却不是素数。

但是经过证明 Miller-Rabin 算法的错误率小于等于 $\dfrac{1}{4}$。若用不同的素数作为 $a$ 测验 $k$ 次，错误率可降低至 $4^{-k}$。

这里还有一个小技巧。如果被测数小于 4759123141 那么只需要测试 3 个底数 2,7,61即可。当然，你测试的越多，范围也就越大。如果你用前7个素数（2,3,5,7,11,13,17）测试，不超过 341550071728320 都是正确的。如果你选用2,3,7,61,24251作为底数，那么 $10^{16}$ 内唯一的伪素数是 46856248255981，是 4840261 的倍数。对于 OIer 来说，已经非常够用了。

Code

当然，对于身为OIer的我们来说，代码肯定是最重要的。

#include<bits/stdc++.h>
#define int long long
using namespace std;
const int Test[5] = {2,3,7,61,24251};
int tsped(int a,int b,int m){
    int ans=0;
    for(;b;b>>=1,(a<<=1)%=m)
        if(b&1)(ans+=a)%=m;
    return ans;
}
int qpow(int a,int b,int m){
    int ans=1;
    for(;b;b>>=1,a=tsped(a,a,m))
        if(b&1)ans=tsped(ans,a,m);
    return ans;
}
bool Miller_Rabin(int n){
    if(n==2)return true;
    if(n<2||!(n&1))return false;
    if(!(n%4840261))return false;
    int m=n-1,q=0;
    while(!(m&1))q++,m>>=1;
    for(int i=0;i<5;i++){
        if(Test[i]==n)continue;
        int a=qpow(Test[i],m,n),b;
        for(int j=1;j<=q;j++){
            b=tsped(a,a,n);
            if(b==1&&a!=1&&a!=n-1)
                return false;
            a=b;
        }
        if(b!=1)return false;
    }
    return true;
}
signed main(){
    int n;
    cin>>n;
    while(n--){
        int a;
        cin>>a;
        cout<<Miller_Rabin(a)<<endl;
    }
    return 0;
}