Miller–Rabin 素性测试

Miller–Rabin 通过幂运算寻找“素数不该有的表现”。它能快速判断 64 位整数是否为素数，不需要先找到因子。

从费马小定理开始

如果 $p$ 是素数， $a$ 不是 $p$ 的倍数，那么

a^{p-1}\equiv1\pmod p.

不满足就能判合数，但满足不一定是素数：

561=3\times11\times17,\qquad 2^{560}\equiv1\pmod{561}.

看看平方是怎么走到 1 的

素数模下， $x^2\equiv1\pmod p$ 可以拆成

(x-1)(x+1)\equiv0\pmod p.

$p$ 是素数，所以它至少整除其中一个因子，也就是 $x\equiv1$ 或 $x\equiv-1$ 。如果一个数既不是 $1$ ，也不是 $-1$ ，平方以后却变成了 $1$ ，模数就一定是合数。

对待测奇数 $n$ ，把指数里的因子 $2$ 全部提出来：

n-1=d\,2^s,\qquad d\text{ 为奇数}.

先算 $a^d\bmod n$ ，再最多平方 $s-1$ 次。素数必须满足：起点为 $1$ ，或途中出现 $n-1$ （即 $-1$ ）。否则判合数。

点选下面三个数，对照一次检查：

多换几个基底

通过一个基底还不能下结论。对 uint64_t，依次检查这七个基底：

\{2,325,9375,28178,450775,9780504,1795265022\}.

在 $0\le n<2^{64}$ 范围内，这组固定基底足以做确定性判定；不能把这个结论直接推广到任意大整数。模乘用 __uint128_t 先接住两个 64 位数的乘积，避免取模前就溢出；这是 GCC / Clang 常见的扩展类型。

代码

展开可直接编译的 C++20 实现

#include <array>
#include <bit>
#include <cstdint>
#include <iostream>

using u64 = std::uint64_t;
using u128 = __uint128_t;

u64 multiplyMod(u64 lhs, u64 rhs, u64 modulus) {
    return static_cast<u64>(static_cast<u128>(lhs) * rhs % modulus);
}

u64 powerMod(u64 base, u64 exponent, u64 modulus) {
    u64 result = 1;
    while (exponent != 0) {
        if (exponent & 1) {
            result = multiplyMod(result, base, modulus);
        }
        base = multiplyMod(base, base, modulus);
        exponent >>= 1;
    }
    return result;
}

bool isPrime(u64 n) {
    if (n < 2) return false;
    if (n % 2 == 0) return n == 2;

    const int s = std::countr_zero(n - 1);
    const u64 d = (n - 1) >> s;
    constexpr std::array<u64, 7> bases{
        2, 325, 9375, 28178, 450775, 9780504, 1795265022
    };

    for (u64 base : bases) {
        base %= n;
        if (base == 0) continue;

        u64 x = powerMod(base, d, n);
        if (x == 1 || x == n - 1) continue;

        bool passed = false;
        for (int r = 1; r < s; ++r) {
            x = multiplyMod(x, x, n);
            if (x == n - 1) {
                passed = true;
                break;
            }
            if (x == 1) return false;
        }
        if (!passed) return false;
    }
    return true;
}

int main() {
    std::ios::sync_with_stdio(false);
    std::cin.tie(nullptr);
    for (u64 n; std::cin >> n;) {
        std::cout << (isPrime(n) ? "Prime" : "Composite") << '\n';
    }
}

每个基底只需要一次快速幂和少量平方；七个基底固定，按一次模乘为常数计，时间是 $O(\log n)$ ，额外空间 $O(1)$ 。

从费马小定理开始

看看平方是怎么走到 1 的

多换几个基底

代码

参考