Miller-Rabin 算法

定义

Miller-Rabin 算法是一种随机化素数测试算法。

有关概念

1. 费马小定理给出了必要条件：当 \(p\) 是质数且 \(p \nmid a\) 时，\(a^{p-1} \equiv 1 \pmod p\)。逆命题不成立，满足该条件的合数称为费马伪素数。
2. 二次探测定理：对素数 \(p\)，同余式 \(x^2 \equiv 1 \pmod p\) 只有两个解，即 \(x \equiv 1 \pmod p\) 或 \(x \equiv -1 \pmod p\)。

求解算法

思路与简单证明

将 \(n-1\) 写成 \(2^s d\)，其中 \(d\) 为奇数。对随机底数 \(a\)，先计算 \(a^d \bmod n\)，再连续平方。若序列既没有出现 \(n-1\)，也不能以 \(1\) 合法结束，则 \(n\) 为合数。对合数而言，单个随机底数误判为素数的概率至多为 \(\frac{1}{4}\)。

实现

下面的实现适用于 64 位无符号整数，使用确定性底数集合。

#include <array>
#include <cstdint>

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

u64 mul_mod(u64 a, u64 b, u64 mod) {
    return static_cast<u128>(a) * b % mod;
}

u64 pow_mod(u64 a, u64 e, u64 mod) {
    u64 result = 1;
    while (e > 0) {
        if (e & 1) result = mul_mod(result, a, mod);
        a = mul_mod(a, a, mod);
        e >>= 1;
    }
    return result;
}

bool is_prime(u64 n) {
    if (n < 2) return false;
    for (u64 p : {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37}) {
        if (n % p == 0) return n == p;
    }
    u64 d = n - 1, s = 0;
    while ((d & 1) == 0) d >>= 1, ++s;
    for (u64 a : {2, 325, 9375, 28178, 450775, 9780504, 1795265022}) {
        if (a % n == 0) continue;
        u64 x = pow_mod(a % n, d, n);
        if (x == 1 || x == n - 1) continue;
        bool composite = true;
        for (u64 r = 1; r < s; ++r) {
            x = mul_mod(x, x, n);
            if (x == n - 1) { composite = false; break; }
        }
        if (composite) return false;
    }
    return true;
}