ホーム > libalgo

確率的素数判定

概要

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

計算量

$n \leq 10^6$ 程度ならエラトステネスのほうが高速

使い方

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

参考文献