可能更好的阅读体验

费马素性测试

费马小定理：
$$a^{p-1} \equiv 1 (\bmod p) [p\in \mathbb{Prime},\gcd(a,p)=1]$$
因为费马小定理的逆命题对于绝大部分的数是正确的，所以我们可以用它来测试素性。
我们可以随机若干个 $a$，对他进行检验：
- 如果 $a^{p-1} \equiv 1(\bmod p)$，那么继续检测。
- 否则其为合数。

练习题目：ARC017A。

int p;
int qpo(int a,int b,const int mod){
    a%=mod;
    int res=1;
    while(b){
        if(b&1) res=(res*a)%mod;
        a=(a*a)%mod;
        b>>=1;
    }
    return res;
}
mt19937 rd(time(NULL));
signed main() {
    in1(p);
    For(i,1,1000000) {
        int a=abs((int)rd());
        if(a%p==0) continue;
        if(qpo(a,p-1,p)!=1) {
            cout<<"NO\n";
            return 0;
        }
    }
    cout<<"YES\n";
    return 0;
}

但是值得注意的是，有这样一类数，如 $561,1105,1729,2465,2821,\cdots$，被称作 Carmichael 数，他们虽然满足对于任意一个 $a$ 都满足 $a^{p-1} \equiv 1 (\bmod p)[\gcd(a,p) = 1]$，但是他们并不是素数，比如 $561=3\times 11 \times 17$，所以我们在随机 $a$ 时可以通过 if(a%p==0) continue; 的方式筛掉大部分的这些数。

当 $\gcd(a,p)=1$ 且 $p$ 为素数时，$a=kp$，所以一定不会出现 $(kp)^{p-1}$。但是 $p$ 不为素数时，$p\mid a$ 不一定代表 $\gcd(a,b)=1$，所以可能出现一个 $a$ 可以筛掉这个数。

Miller-Rabin 素性测试

https://oi-wiki.org/math/number-theory/prime/#millerrabin-%E7%B4%A0%E6%80%A7%E6%B5%8B%E8%AF%95

https://www.luogu.com.cn/article/jr27vtto

二次探测定理：$x^2\equiv 1(\bmod p)$ 且 $p\ge 3$ 且 $p$ 为素数时，$x$ 的解为：$x\equiv 1(\bmod p)$ 或 $x\equiv p-1(\mod p)$。

$$\because x^2 \equiv 1(\bmod p)$$
$$\therefore x^2-1 \equiv 0(\bmod p)$$
$$\therefore (x-1)(x+1) \equiv 0 (\bmod p)$$
$$x-1 \mid p 或 x+1 \mid p$$
$$x_1\equiv 1(\bmod p),x_2\equiv p-1(\bmod p)$$

设 $p-1=2^tu$，那么根据费马小定理，我们有 $a^{p-1}\equiv (a^u)^{2^t} (\bmod p)$，用快速幂计算 $a^u$。然后重复执行 $t$ 次平方。因为每一次平方都可以用二次探测定理检验，遂最多计算 $t$ 次平方。

通过选取一些固定的数字，也可以实现一定范围内的正确性。
在考场上选取 $\{2,3,5,7,11,13,17,19,23,29,31,37\}$ 可以判断 $2^{78}$ 以内的数字。
如果你想知道前 $n$ 个素数可以判断多少数字，参见A014233。

练习题目：PON - Prime or Not。

#define int __int128
int t,p;
int tst[15]={0,2,3,5,7,11,13,17,19,23,29,31,37};
int qpo(int a,int b,int mod) {
    int res=1;
    for(;b;b>>=1,a=(a*a)%mod) if(b&1) res=(res*a)%mod;
    return res;
}
void divide(int p,int &x,int &u) {
    while(p%(1<<(x+1))==0) x++;
    u=p/(1<<x);
}
bool check(int p) {
    if(p<=40) {
        For(i,1,12) if(p==tst[i]) return 1;
        return 0;
    }
    int t=0,u=0;
    divide(p-1,t,u);
    For(i,1,12) {
        int a=tst[i];
        a=qpo(a,u,p);
        if(a==1) continue;
        bool flg=0;
        For(j,1,t) {
            if(a==p-1) {
                flg=1;
                break;
            }
            a=(a*a)%p;
        }
        if(!flg) return 0;
    }
    return 1;
}