素数判定 (Math/primality_test.hpp)

View this file on GitHub
View document part on GitHub
Last update: 2025-04-13 19:22:19+09:00
Include: #include "Math/primality_test.hpp"

PrimalityTest

Miller-Rabin素数判定法を用いた素数判定を行います。
内部での累乗の計算を128ビット整数を用いた繰り返し二乗法で行っているため、やや速度が遅くなっています。
モンゴメリ乗算をライブラリに追加したら累乗計算の部分は置き換えたい。

PrimalityTest::is_prime

static bool PrimalityTest::is_prime(long long N)

整数Nが素数であるかどうかを判定します。

制約

$0 \leq N$

計算量

$O(\log N)$

Depends on

Required by

高速素因数分解 (Math/fast_factorizer.hpp)

Verified with

Code

#include <bits/stdc++.h>
using namespace std;

#include "Math/mod_power.hpp"

struct PrimalityTest
{
private:
    static bool miller_rabin(long long N, const vector<long long> &A)
    {
        const int s = countr_zero((unsigned long long)N - 1);
        const long long d = (N - 1) >> s;

        auto is_composite = [&](long long a) -> bool
        {
            // Fermat の小定理が成り立たない
            if (a % N == 0)
            {
                return false;
            }

            long long x = mod_pow<__int128_t>(a, d, N);
            if (x == 1 || x == N - 1)
            {
                return false;
            }

            for (int t = 0; t < s - 1; ++t)
            {
                x = __int128_t(x) * x % N;
                if (x == N - 1)
                {
                    return false;
                }
            }

            return true;
        };

        return none_of(A.begin(), A.end(), is_composite);
    }

public:
    static bool is_prime(long long N)
    {
        assert(0 <= N);
        if (N <= 2)
        {
            return N == 2;
        }
        if (N % 2 == 0)
        {
            return false;
        }

        return miller_rabin(N, {2, 325, 9375, 28178, 450775, 9780504, 1795265022});
    }
};

#line 1 "Math/primality_test.hpp"
#include <bits/stdc++.h>
using namespace std;

#line 2 "Math/mod_power.hpp"
using namespace std;

#line 2 "Other/fast_power.hpp"
using namespace std;

template <class S>
S fast_pow(S x, long long n, function<S(S, S)> mul, function<S()> e)
{
    assert(0 <= n);

    S ans = e();

    while (n)
    {
        if (n & 1)
        {
            ans = mul(ans, x);
        }
        x = mul(x, x);
        n >>= 1;
    }

    return ans;
}
#line 5 "Math/mod_power.hpp"

template <typename T>
enable_if_t<is_integral_v<T> || is_same_v<T, __int128_t>, T>
mod_pow(T x, T n, T mod)
{
    assert(0 <= n);
    assert(0 < mod);
    assert(mod <= numeric_limits<T>::max() / mod);

    x %= mod;
    if (x < 0)
    {
        x += mod;
    }

    auto mul = [&](T a, T b) -> T
    {
        return (a * b) % mod;
    };

    auto e = [&]() -> T
    {
        return 1;
    };

    return fast_pow<T>(x, n, mul, e);
}
#line 5 "Math/primality_test.hpp"

struct PrimalityTest
{
private:
    static bool miller_rabin(long long N, const vector<long long> &A)
    {
        const int s = countr_zero((unsigned long long)N - 1);
        const long long d = (N - 1) >> s;

        auto is_composite = [&](long long a) -> bool
        {
            // Fermat の小定理が成り立たない
            if (a % N == 0)
            {
                return false;
            }

            long long x = mod_pow<__int128_t>(a, d, N);
            if (x == 1 || x == N - 1)
            {
                return false;
            }

            for (int t = 0; t < s - 1; ++t)
            {
                x = __int128_t(x) * x % N;
                if (x == N - 1)
                {
                    return false;
                }
            }

            return true;
        };

        return none_of(A.begin(), A.end(), is_composite);
    }

public:
    static bool is_prime(long long N)
    {
        assert(0 <= N);
        if (N <= 2)
        {
            return N == 2;
        }
        if (N % 2 == 0)
        {
            return false;
        }

        return miller_rabin(N, {2, 325, 9375, 28178, 450775, 9780504, 1795265022});
    }
};