miller-rabin primality test
(number/miller-rabin.hpp)

• View this file on GitHub
• Last update: 2021-08-30 04:35:37+09:00
• Include: #include "number/miller-rabin.hpp"

Depends on

• modulus/modpow.hpp
• utils/macros.hpp

Code

#pragma once
#include <cassert>
#include <random>
#include "../utils/macros.hpp"
#include "../modulus/modpow.hpp"

/**
 * @brief miller-rabin primality test
 * @note $O(k \log n)$
 */
template <class RandomEngine>
bool is_prime(int64_t n, int iteration, RandomEngine& gen) {
    assert (0 <= n);
    if (n == 2) return true;
    if (n == 1 or n % 2 == 0) return false;
    const int64_t d = (n - 1) >> __builtin_ctzll(n - 1);  // remove trailing zeros
    std::uniform_int_distribution<int64_t> dist(1, n - 2);  // [l, r]
    REP (dummy, iteration) {
        int64_t a = dist(gen);
        int64_t t = d;
        int64_t y = modpow(a, t, n);
        while (t != n - 1 and y != 1 and y != n - 1) {
            y = y * y % n;
            t *= 2;
        }
        if (y != n - 1 and t % 2 == 0) return false;
    }
    return true;
}

bool is_prime(int64_t n) {
    static std::default_random_engine gen = std::default_random_engine(std::random_device()());
    return is_prime(n, 20, gen);
}

#line 2 "number/miller-rabin.hpp"
#include <cassert>
#include <random>
#line 2 "utils/macros.hpp"
#define REP(i, n) for (int i = 0; (i) < (int)(n); ++ (i))
#define REP3(i, m, n) for (int i = (m); (i) < (int)(n); ++ (i))
#define REP_R(i, n) for (int i = (int)(n) - 1; (i) >= 0; -- (i))
#define REP3R(i, m, n) for (int i = (int)(n) - 1; (i) >= (int)(m); -- (i))
#define ALL(x) std::begin(x), std::end(x)
#line 3 "modulus/modpow.hpp"
#include <cstdint>

inline int32_t modpow(uint_fast64_t x, uint64_t k, int32_t MOD) {
    assert (/* 0 <= x and */ x < (uint_fast64_t)MOD);
    uint_fast64_t y = 1;
    for (; k; k >>= 1) {
        if (k & 1) (y *= x) %= MOD;
        (x *= x) %= MOD;
    }
    assert (/* 0 <= y and */ y < (uint_fast64_t)MOD);
    return y;
}
#line 6 "number/miller-rabin.hpp"

/**
 * @brief miller-rabin primality test
 * @note $O(k \log n)$
 */
template <class RandomEngine>
bool is_prime(int64_t n, int iteration, RandomEngine& gen) {
    assert (0 <= n);
    if (n == 2) return true;
    if (n == 1 or n % 2 == 0) return false;
    const int64_t d = (n - 1) >> __builtin_ctzll(n - 1);  // remove trailing zeros
    std::uniform_int_distribution<int64_t> dist(1, n - 2);  // [l, r]
    REP (dummy, iteration) {
        int64_t a = dist(gen);
        int64_t t = d;
        int64_t y = modpow(a, t, n);
        while (t != n - 1 and y != 1 and y != n - 1) {
            y = y * y % n;
            t *= 2;
        }
        if (y != n - 1 and t % 2 == 0) return false;
    }
    return true;
}

bool is_prime(int64_t n) {
    static std::default_random_engine gen = std::default_random_engine(std::random_device()());
    return is_prime(n, 20, gen);
}

Back to top page