cplib/miller_rabin.hpp

View this file on GitHub
Last update: 2025-12-24 16:17:55+09:00
Include: #include "cplib/miller_rabin.hpp"

Depends on

cplib/pow_mod_ull.hpp

Required by

cplib/is_prime.hpp

Verified with

Code

#pragma once

#include <bit>
#include <cstdint>
#include <ranges>
#include <vector>
#include "cplib/pow_mod_ull.hpp"

namespace cplib {

// n が素数であるかをミラーラビン法で判定する
// vector で渡した整数を証人とする
bool miller_rabin(const uint64_t n, const std::vector<uint64_t>& wits) {
    if (n < 2) return false;
    if (n == 2) return true;
    if (n % 2 == 0) return false;

    const int s = std::countr_zero(n - 1);
    const uint64_t d = (n - 1) >> s;

    for (const uint64_t a :
         wits | std::views::transform([n](uint64_t x) { return x % n; })) {
        if (a == 0) continue;

        uint64_t tmp = pow_mod_ull(a, d, n);
        if (tmp == 1) continue;

        for (int i = 0; tmp != n - 1; ++i) {
            if (i == s - 1) return false;
            tmp = pow_mod_ull(tmp, 2, n);
        }
    }

    return true;
}

}  // namespace cplib

#line 2 "cplib/miller_rabin.hpp"

#include <bit>
#include <cstdint>
#include <ranges>
#include <vector>
#line 2 "cplib/pow_mod_ull.hpp"

#line 4 "cplib/pow_mod_ull.hpp"

namespace cplib {

// (x ^ n) mod m を繰り返し二乗法で求める
// x >= 0, n >= 0, m >= 1
uint64_t pow_mod_ull(uint64_t x, uint64_t n, const uint64_t m) {
    uint64_t result = ((n % 2 == 1) ? x : 1) % m;
    for (n /= 2; n > 0; n /= 2) {
        x = static_cast<uint64_t>(__uint128_t(x) * __uint128_t(x) %
                                  __uint128_t(m));
        if (n % 2 == 1) {
            result = static_cast<uint64_t>(__uint128_t(result) *
                                           __uint128_t(x) % __uint128_t(m));
        }
    }
    return result;
}

}  // namespace cplib
#line 8 "cplib/miller_rabin.hpp"

namespace cplib {

// n が素数であるかをミラーラビン法で判定する
// vector で渡した整数を証人とする
bool miller_rabin(const uint64_t n, const std::vector<uint64_t>& wits) {
    if (n < 2) return false;
    if (n == 2) return true;
    if (n % 2 == 0) return false;

    const int s = std::countr_zero(n - 1);
    const uint64_t d = (n - 1) >> s;

    for (const uint64_t a :
         wits | std::views::transform([n](uint64_t x) { return x % n; })) {
        if (a == 0) continue;

        uint64_t tmp = pow_mod_ull(a, d, n);
        if (tmp == 1) continue;

        for (int i = 0; tmp != n - 1; ++i) {
            if (i == s - 1) return false;
            tmp = pow_mod_ull(tmp, 2, n);
        }
    }

    return true;
}

}  // namespace cplib