確率的素数判定

概要

$p$ が素数である必要条件を何種類かの候補 $a$ でテストして，$p$ が素数だと確信を得ていくアルゴリズムの一つにミラー・ラビンテスト (Miller-Rabin Test) がある．テストに失敗する $a$ のことを $p$ が合成数であることの証人 (witness) という．普通，証人は $[2,n-1]$ から何度か乱択するが，$p \leq 2^{32}, 2^{64}$ の場合に限ると厳密な判定ができる小さな証人の集合が存在する．

計算量

$O(k \log(n)^3)$．ただし $k$ は異なる $a$ の個数．実験はしていないが，定数が大きいのでエラトステネスの篩や試し割りの方が高速なことも多いと思う．

使い方

is_prime(n): $n$ が素数か判定する．$n < 2^{63}$ ならば厳密な判定を行う．__uint128_t が使えないときは modmul をコメントアウトされているビット演算版に切り替える．

実装

using u32 = unsigned int;
using u64 = unsigned long long;
using u128 = __uint128_t;

template <class Uint, class BinOp>
bool is_prime_impl(const Uint &n, const Uint *witness, BinOp modmul) {
    if (n == 2) return true;
    if (n < 2 || n % 2 == 0) return false;
    const Uint m = n - 1, d = m / (m & -m);
    auto modpow = [&](Uint a, Uint b) {
        Uint res = 1;
        for (; b; b /= 2) {
            if (b & 1) res = modmul(res, a);
            a = modmul(a, a);
        }
        return res;
    };
    auto suspect = [&](Uint a, Uint t) {
        a = modpow(a, t);
        while (t != n - 1 && a != 1 && a != n - 1) {
            a = modmul(a, a);
            t = modmul(t, 2);
        }
        return a == n - 1 || t % 2 == 1;
    };
    for (const Uint *w = witness; *w; w++) {
        if (*w % n != 0 && !suspect(*w, d)) return false;
    }
    return true;
}

bool is_prime(const u128 &n) {
    assert(n < 1ULL << 63);
    if (n < 1ULL << 32) {
        // n < 2^32
        constexpr u64 witness[] = {2, 7, 61, 0};
        auto modmul = [&](u64 a, u64 b) { return a * b % n; };
        return is_prime_impl<u64>(n, witness, modmul);
    } else {
        // n < 2^63
        constexpr u128 witness[] = {2, 325, 9375, 28178, 450775, 9780504, 1795265022, 0};
        // if u128 is available
        auto modmul = [&](u128 a, u128 b) { return a * b % n; };
        // otherwise
        // auto modmul = [&](u64 a, u64 b) {
        //     u64 res = 0;
        //     for (; b; b /= 2) {
        //         if (b & 1) res = (res + a) % n;
        //         a = (a + a) % n;
        //     }
        //     return res;
        // };
        return is_prime_impl<u128>(n, witness, modmul);
    }
}

検証

AOJ ALDS11C http://judge.u-aizu.ac.jp/onlinejudge/review.jsp?rid=2169070

SPOJ PON http://ideone.com/IkXM3S

参考文献

ミラー–ラビン素数判定法 - Wikipedia https://ja.wikipedia.org/wiki/%E3%83%9F%E3%83%A9%E3%83%BC%E2%80%93%E3%83%A9%E3%83%93%E3%83%B3%E7%B4%A0%E6%95%B0%E5%88%A4%E5%AE%9A%E6%B3%95
Fast Primality Testing for Integers That Fit into a Machine Word http://ceur-ws.org/Vol-1326/020-Forisek.pdf