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;
}