Library/Math/MillerRabin.hpp

• View this file on GitHub
• Last update: 2026-06-01 13:41:50+09:00
• Include: #include "Library/Math/MillerRabin.hpp"

Depends on

• Library/Common.hpp
• Library/Math/Montgomery.hpp

Verified with

• verify/LC-PrimalityTest.test.cpp

Code

#include "Montgomery.hpp"

bool MillerRabin(const long long N){
    if(N <= 1) return false;
    if(N == 2) return true;
    if(!(N & 1)) return false;
    vector<long long> as;
    if(N < 4759123141LL) as = {2, 7, 61};
    else as = {2, 325, 9375, 28178, 450775, 9780504, 1795265022};
    Montgomery::SetMod(N);
    long long s = 0, d = N - 1;
    while(!(d & 1)) ++s, d >>= 1;
    for(const auto &a : as){
        Montgomery x = Montgomery(a).Power(d);
        if(x != 1 && x != 0){
            long long t;
            for(t = 0; t < s; ++t){
                if(x == N - 1) break;
                x *= x;
            }
            if(t == s) return false;
        }
    }
    return true;
}

#line 2 "Library/Common.hpp"

/**
 * @file Common.hpp
 */

#include <algorithm>
#include <array>
#include <bit>
#include <bitset>
#include <cassert>
#include <cmath>
#include <cstdint>
#include <deque>
#include <functional>
#include <iomanip>
#include <iostream>
#include <limits>
#include <map>
#include <numeric>
#include <queue>
#include <set>
#include <stack>
#include <string>
#include <tuple>
#include <utility>
#include <vector>
using namespace std;

using ll = int64_t;
using ull = uint64_t;

constexpr const ll INF = (1LL << 62) - (3LL << 30) - 1;
#line 2 "Library/Math/Montgomery.hpp"

struct Montgomery{
    using mont = Montgomery;
    using u64 = uint64_t;
    using u128 = __uint128_t;

    static void SetMod(const u64 m){
        assert(m < (1LL << 62));
        assert(m & 1);
        mod = m;
        inv_mod = m;
        for(int i = 0; i < 5; ++i){
            inv_mod *= 2 - m * inv_mod;
        }
        r2 = -u128(m) % m;
    }

    u64 val;

    Montgomery() : val(0){}

    Montgomery(long long x) : val(Reduct((u128(x) + mod) * r2)){}

    inline u64 Get() const {
        u64 res = Reduct(val);
        return res >= mod ? res - mod : res;
    }

    mont &operator+=(const mont &r){
        if((val += r.val) >= 2 * mod) val -= 2 * mod;
        return (*this);
    }

    mont &operator-=(const mont &r){
        if((val += 2 * mod - r.val) >= 2 * mod) val -= 2 * mod;
        return (*this);
    }

    mont &operator*=(const mont &r){
        val = Reduct(u128(val) * r.val);
        return (*this);
    }

    mont &operator/=(const mont &r){
        (*this) *= r.Inverse();
        return (*this);
    }

    mont Inverse() const {
        return Power(mod - 2);
    }

    mont Power(long long k) const {
        mont ret(1), mul(*this);
        while(k){
            if(k & 1) ret *= mul;
            mul *= mul;
            k >>= 1;
        }
        return ret;
    }

    mont operator-() const {
        return mont() - mont(*this);
    }

    mont operator+(const mont &r) const {
        return mont(*this) += r;
    }

    mont operator-(const mont &r) const {
        return mont(*this) -= r;
    }

    mont operator*(const mont &r) const {
        return mont(*this) *= r;
    }

    bool operator==(const mont &r) const {
        return (val >= mod ? val - mod : val) == (r.val >= mod ? r.val - mod : r.val);
    }

    bool operator!=(const mont &r) const {
        return (val >= mod ? val - mod : val) != (r.val >= mod ? r.val - mod : r.val);
    }

    private:
    static u64 mod;
    static u64 inv_mod;
    static u64 r2;

    static u64 Reduct(const u128 x){
        return (x + u128(u64(x) * u64(-inv_mod)) * mod) >> 64;
    }
};

typename Montgomery::u64 Montgomery::mod, Montgomery::inv_mod, Montgomery::r2;
#line 2 "Library/Math/MillerRabin.hpp"

bool MillerRabin(const long long N){
    if(N <= 1) return false;
    if(N == 2) return true;
    if(!(N & 1)) return false;
    vector<long long> as;
    if(N < 4759123141LL) as = {2, 7, 61};
    else as = {2, 325, 9375, 28178, 450775, 9780504, 1795265022};
    Montgomery::SetMod(N);
    long long s = 0, d = N - 1;
    while(!(d & 1)) ++s, d >>= 1;
    for(const auto &a : as){
        Montgomery x = Montgomery(a).Power(d);
        if(x != 1 && x != 0){
            long long t;
            for(t = 0; t < s; ++t){
                if(x == N - 1) break;
                x *= x;
            }
            if(t == s) return false;
        }
    }
    return true;
}

Back to top page