1. 是什么

用来在 $O(\log p)$的复杂度内检测 $p$是否为素数

是一种随机化玄学算法，不一定保证正确

正确率后面会讲

2. 算法

首先由费马小定理可知：

若 $p$为素数，那么对于任意一个 $a \in [1, p – 1]$都有：

$$a ^ {p – 1} \equiv 1 \mod p$$

因此可以任取一个 $a \in [1, p – 1]$，快速幂检测是否满足该式，若不满足则一定为合数，若满足并不一定是质数，你可能取少了，因此要多取几组

但是现在有一个问题，有的比较奇妙的 $p$是合数，但是对于任意的 $a \in [1, p – 1]$，都满足 $a ^ {p – 1} \equiv 1 \mod p$，那怎么办呢？

于是我们引入 “二次探测定理”：

如果 p 是奇素数，则 $x ^ 2 \equiv 1 \mod p$的解为 $x = 1$ 或 $x = p − 1$（指的是在 $\mod p$意义下的解）

因为 $a ^ {p – 1} \equiv 1 \mod p$，若要满足二次探测定理，则需要满足 $a ^ \frac {p – 1} 2 \equiv 1 \mod p$或 $\equiv p – 1 \mod p$

因此当 $p – 1$为偶数时，令 $i = \frac {p – 1} 2$，若 $a ^ i \not \equiv 1 \mod p$且 $a ^ i \not \equiv p – 1 \mod p$则说明 $p$不是素数

若 $a ^ i \equiv 1 \mod p$则继续判断，即判断 $i$是否为偶数，如果为偶数则判断 $a ^ \frac i 2 \mod p$的值是否等于 $1$或 $p – 1$…

直到 $i$为奇数（或者 $a ^ i \equiv p – 1 \mod p$），如果此时依然没有出现过二次探测失败的情况则说明二次探测成功

这样可以使得 Miller-Rabin 的正确率达到 $1 – \frac 1 {4 ^ k}$，其中 $k$表示你取了多少个 $a \in [1, p – 1]$

（但令我吃惊的是只写了费马小定理的 Miller-Rabin 居然在 hihoCoder 上能过）

（当然在 LOJ 上是过不了的）

你可以算一下，当 $k = 10$时 Miller-Rabin 的正确率为 $0.999999046326$，要是这都能 WA。。。可以去买⑥合采了

如果直接地先跑费马小定理，再二次探测，每次探测都跑一下快速幂，那么复杂度就是 $O(\log ^ 2 n)$的（带个 $k$的常数）

如果你想达到 $O(\log n)$的复杂度，就先算出一个 $k$和一个 $t$，使得 $2 ^ t k = p – 1$且 $2 \not \mid k$，然后快速幂算出 $a ^ k$，令 $s = a ^ k$，接着循环 $k$次，每次判断 $s ^ 2 \mod p$是否等于 $1$，若 $s ^ 2 \mod p$等于 $1$则判断 $s$是否等于 $1$或者 $p – 1$，若不等于则二次探测失败，否则令 $s = s ^ 2$，进行下一轮判断

完成 $t$轮二次探测后 $s = a ^ {p – 1}$，再进行费马小定理就行了

复杂度 $O(\log n)$

如果数据较大（$n \leq 10 ^ {18}$），则需要快速乘，就是多个 $\log$而已，如果不想多这个 $\log$就用 $O(1)$快速乘（参见下方代码里的 mul 函数）

例题

LOJ – 143 质数判定

代码：

#include <bits/stdc++.h>

using namespace std;

typedef long long LL;

LL mul(LL a, LL b, LL M)
{
    a = a % M, b = b % M;
    return ((a * b - (LL)(((long double) a * b + 0.5) / M) * M) % M + M) % M;
}

LL n;

LL qpow(LL a, LL b, LL c)
{
    LL res = 1;
    while (b)
    {
        if (b & 1) res = mul(res, a, c);
        a = mul(a, a, c), b >>= 1;
    }
    return res % c;
}

bool mr(LL p)
{
    if (p < 2) return 0;
    if (p == 2) return 1;
    if (!(p & 1)) return 0;
    LL k = p - 1; int t = 0;
    while (!(k & 1)) k >>= 1, t++;
    for (int C = 1; C <= 10; C += 1)
    {
        LL a = rand() % (p - 2) + 2, b = qpow(a, k, p), c = b;
        for (int i = 1; i <= t; i += 1)
        {
            b = mul(b, b, p);
            if (b == 1 && c != 1 && c != p - 1) return 0;
            c = b;
        }
        if (b != 1) return 0;
    }
    return 1;
}

int main(int argc, char const* argv[])
{
    srand(5128957);
    while (~scanf("%lld", &n))
        if (mr(n)) puts("Y");
        else puts("N");
    return 0;
}