#865710 (C++23) No.8030 ミラー・ラビン素数判定法のテスト

提出ソース
結果

問題	No.8030 ミラー・ラビン素数判定法のテスト
ユーザー	AC2K
提出日時	2023-05-03 22:56:40
言語	C++23 (gcc 13.3.0 + boost 1.87.0)
結果	AC
実行時間	44 ms / 9,973 ms
コード長	11,852 bytes
コンパイル時間	764 ms
コンパイル使用メモリ	83,736 KB
実行使用メモリ	5,248 KB
最終ジャッジ日時	2024-11-21 22:14:28
合計ジャッジ時間	1,507 ms
ジャッジサーバーID （参考情報）	judge4 / judge2
このコードへのチャレンジ
（要ログイン）
ファイルパターン	結果
other	AC * 10
権限があれば一括ダウンロードができます
ソースコード

raw source code
#line 1 "test/yuki/No3030.test.cpp"
#include<iostream>
#line 2 "src/math/dynamic_modint.hpp"
#include <cassert>
#line 2 "src/internal/barrett.hpp"
namespace kyopro {
namespace internal {
/// @brief barrett reduction
/// @ref https://github.com/atcoder/ac-library/blob/master/atcoder/internal_math.hpp
class barrett {
    using u32 = uint32_t;
    using u64 = uint64_t;

    u64 m;
    u64 im;

public:
    explicit barrett() = default;
    explicit barrett(u64 m_)
        : m(m_), im((u64)(long double)static_cast<u64>(-1) / m_ + 1) {}

    u64 get_mod() const { return m; }
    constexpr u64 reduce(int64_t a) const {
        if (a < 0) return m - reduce(-a);
        u64 q = ((__uint128_t)a * im) >> 64;
        a -= m * q;
        if (a >= m) a -= m;
        return a;
    }
    constexpr u64 mul(u64 a, u64 b) const {
        if (a == 0 || b == 0) {
            return 0;
        }
        u64 z = a;
        z *= b;
        u64 x = (u64)(((__uint128_t)z * im) >> 64);

        u32 v = (u32)(z - x * m);

        if (v >= m) v += m;
        return v;
    }
};
};  // namespace internal
};  // namespace kyopro
#line 3 "src/internal/montgomery.hpp"
#include <limits>
#include <numeric>
#line 5 "src/internal/type_traits.hpp"
#include <typeinfo>
namespace kyopro {
namespace internal {
/// @ref https://qiita.com/kazatsuyu/items/f8c3b304e7f8b35263d8
template <typename... Args> struct first_enabled {};

template <typename T, typename... Args>
struct first_enabled<std::enable_if<true, T>, Args...> {
    using type = T;
};
template <typename T, typename... Args>
struct first_enabled<std::enable_if<false, T>, Args...>
    : first_enabled<Args...> {};
template <typename T, typename... Args> struct first_enabled<T, Args...> {
    using type = T;
};

template <typename... Args>
using first_enabled_t = typename first_enabled<Args...>::type;

template <int dgt> struct int_least {
    static_assert(dgt <= 128);
    using type = first_enabled_t<std::enable_if<dgt <= 8, __int8_t>,
                                 std::enable_if<dgt <= 16, __int16_t>,
                                 std::enable_if<dgt <= 32, __int32_t>,
                                 std::enable_if<dgt <= 64, __int64_t>,
                                 std::enable_if<dgt <= 128, __int128_t> >;
};
template <int dgt> struct uint_least {
    static_assert(dgt <= 128);
    using type = first_enabled_t<std::enable_if<dgt <= 8, __uint8_t>,
                                 std::enable_if<dgt <= 16, __uint16_t>,
                                 std::enable_if<dgt <= 32, __uint32_t>,
                                 std::enable_if<dgt <= 64, __uint64_t>,
                                 std::enable_if<dgt <= 128, __uint128_t> >;
};

template <int dgt> using int_least_t = typename int_least<dgt>::type;
template <int dgt> using uint_least_t = typename uint_least<dgt>::type;

template <typename T>
using double_size_uint_t = uint_least_t<2 * std::numeric_limits<T>::digits>;

template <typename T>
using double_size_int_t = int_least_t<2 * std::numeric_limits<T>::digits>;
};  // namespace internal
};  // namespace kyopro
#line 6 "src/internal/montgomery.hpp"
namespace kyopro {
namespace internal {
using u32 = uint32_t;
using u64 = uint64_t;
using i32 = int32_t;
using i64 = int64_t;
using u128 = __uint128_t;
using i128 = __int128_t;
/// @brief MontgomeryReduction
/// @ref 
template <typename T>
class Montgomery {
    static constexpr int lg = std::numeric_limits<T>::digits;
    using LargeT = internal::double_size_uint_t<T>;
    T mod, r, r2, minv;
    T inv() {
        T t = 0, res = 0;
        for (int i = 0; i < lg; ++i) {
            if (~t & 1) {
                t += mod;
                res += static_cast<T>(1) << i;
            }
            t >>= 1;
        }
        return res;
    }

public:
    Montgomery() = default;
    constexpr T get_mod() { return mod; }
    constexpr int get_lg() { return lg; }

    void set_mod(T m) {
        assert(m > 0);
        assert(m & 1);

        mod = m;

        r = (-static_cast<T>(mod)) % mod;
        r2 = (-static_cast<LargeT>(mod)) % mod;
        minv = inv();
    }

    T reduce(LargeT x) const {
        u64 res =
            (x + static_cast<LargeT>(static_cast<T>(x) * minv) * mod) >> lg;

        if (res >= mod) res -= mod;
        return res;
    }

    T generate(LargeT x) { return reduce(x * r2); }

    T mult(T x, T y) { return reduce((LargeT)x * y); }
};
};  // namespace internal
};  // namespace kyopro
#line 6 "src/math/dynamic_modint.hpp"
namespace kyopro {
/// @note mod は32bitじゃないとバグる
template <int id = -1>
class barrett_modint {
    using u32 = uint32_t;
    using u64 = uint64_t;

    using i32 = int32_t;
    using i64 = int64_t;
    using br = internal::barrett;

    static br brt;
    static u32 mod;
    u32 v;  // value
public:
    static void set_mod(u32 mod_) {
        brt = br(mod_);
        mod = mod_;
    }

public:
    explicit constexpr barrett_modint() : v(0) {
        assert(mod);
    }
    explicit constexpr barrett_modint(i64 v_) : v(brt.reduce(v_)) {
        assert(mod);
    }

    u32 val() const { return v; }
    static u32 get_mod() { return mod; }
    using mint = barrett_modint<id>;

    constexpr mint& operator=(i64 r) {
        v = brt.reduce(r);
        return (*this);
    }
    constexpr mint& operator+=(const mint& r) {
        v += r.v;
        if (v >= mod) {
            v -= mod;
        }
        return (*this);
    }
    constexpr mint& operator-=(const mint& r) {
        v += mod - r.v;
        if (v >= mod) {
            v -= mod;
        }

        return (*this);
    }
    constexpr mint& operator*=(const mint& r) {
        v = brt.mul(v, r.v);
        return (*this);
    }

    constexpr mint operator+(const mint& r) const { return mint(*this) += r; }
    constexpr mint operator-(const mint& r) const { return mint(*this) -= r; }
    constexpr mint operator*(const mint& r) const { return mint(*this) *= r; }

    constexpr mint& operator+=(i64 r) { return (*this) += mint(r); }
    constexpr mint& operator-=(i64 r) { return (*this) -= mint(r); }
    constexpr mint& operator*=(i64 r) { return (*this) *= mint(r); }

    friend mint operator+(i64 l, const mint& r) { return mint(l) += r; }
    friend mint operator+(const mint& l, i64 r) { return mint(l) += r; }
    friend mint operator-(i64 l, const mint& r) { return mint(l) -= r; }
    friend mint operator-(const mint& l, i64 r) { return mint(l) -= r; }
    friend mint operator*(i64 l, const mint& r) { return mint(l) *= r; }
    friend mint operator*(const mint& l, i64 r) { return mint(l) += r; }

    friend std::ostream& operator<<(std::ostream& os, const mint& mt) {
        os << mt.val();
        return os;
    }
    friend std::istream& operator>>(std::istream& is, mint& mt) {
        i64 v_;
        is >> v_;
        mt = v_;
        return is;
    }
    template <typename T>
    mint pow(T e) const {
        mint res(1), base(*this);

        while (e) {
            if (e & 1) {
                res *= base;
            }
            e >>= 1;
            base *= base;
        }
        return res;
    }
    mint inv() const { return pow(mod - 2); }

    mint& operator/=(const mint& r) { return (*this) *= r.inv(); }
    mint operator/(const mint& r) const { return mint(*this) *= r.inv(); }
    mint& operator/=(i64 r) { return (*this) /= mint(r); }
    friend mint operator/(const mint& l, i64 r) { return mint(l) /= r; }
    friend mint operator/(i64 l, const mint& r) { return mint(l) /= r; }
};
};  // namespace kyopro
template <int id>
typename kyopro::barrett_modint<id>::u32 kyopro::barrett_modint<id>::mod;
template <int id>
typename kyopro::barrett_modint<id>::br kyopro::barrett_modint<id>::brt;

namespace kyopro {
template <typename T, int id = -1>
class dynamic_modint {
    using LargeT = internal::double_size_uint_t<T>;
    static T mod;
    static internal::Montgomery<T> mr;

public:
    static void set_mod(T mod_) {
        mr.set_mod(mod_);
        mod = mod_;
    }

    static T get_mod() { return mod; }

private:
    T v;

public:
    dynamic_modint(T v_ = 0) {
        assert(mod);
        v = mr.generate(v_);
    }
    T val() const { return mr.reduce(v); }

    using mint = dynamic_modint<T, id>;
    mint& operator+=(const mint& r) {
        v += r.v;
        if (v >= mr.get_mod()) {
            v -= mr.get_mod();
        }

        return (*this);
    }

    mint& operator-=(const mint& r) {
        v += mr.get_mod() - r.v;
        if (v >= mr.get_mod) {
            v -= mr.get_mod();
        }

        return (*this);
    }

    mint& operator*=(const mint& r) {
        v = mr.mult(v, r.v);
        return (*this);
    }

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

    mint& operator=(const T& v_) {
        (*this) = mint(v_);
        return (*this);
    }

    friend std::ostream& operator<<(std::ostream& os, const mint& mt) {
        os << mt.val();
        return os;
    }
    friend std::istream& operator>>(std::istream& is, mint& mt) {
        T v_;
        is >> v_;
        mt = v_;
        return is;
    }
    template <typename P>
    mint pow(P e) const {
        assert(e >= 0);
        mint res(1), base(*this);

        while (e) {
            if (e & 1) {
                res *= base;
            }
            e >>= 1;
            base *= base;
        }
        return res;
    }
    mint inv() const { return pow(mod - 2); }

    mint& operator/=(const mint& r) { return (*this) *= r.inv(); }
    mint operator/(const mint& r) const { return mint(*this) *= r.inv(); }
    mint& operator/=(T r) { return (*this) /= mint(r); }
    friend mint operator/(const mint& l, T r) { return mint(l) /= r; }
    friend mint operator/(T l, const mint& r) { return mint(l) /= r; }
};
};  // namespace kyopro
template <typename T, int id>
T kyopro::dynamic_modint<T, id>::mod;
template <typename T, int id>
kyopro::internal::Montgomery<T> kyopro::dynamic_modint<T, id>::mr;

/// @brief dynamic modint
/// @docs docs/math/dynamic_modint.md
#line 3 "src/math/miller.hpp"
namespace kyopro {
namespace miller {
using i128 = __int128_t;
using u128 = __uint128_t;
using u64 = uint64_t;
using u32 = uint32_t;

using i128 = __int128_t;
using u128 = __uint128_t;
using u64 = uint64_t;
using u32 = uint32_t;

template <typename T, typename mint,const int bases[],int length>
constexpr bool miller_rabin(T n) {
    T d = n - 1;

    while (~d & 1) {
        d >>= 1;
    }

    const T rev = n - 1;
    if (mint::get_mod() != n) {
        mint::set_mod(n);
    }
    for (int i = 0; i < length; ++i) {
        if (n <= bases[i]) {
            return true;
        }
        T t = d;
        mint y = mint(bases[i]).pow(t);
        while (t != n - 1 && y.val() != 1 && y.val() != rev) {
            y *= y;
            t <<= 1;
        }

        if (y.val() != rev && (~t & 1)) return false;
    }
    return true;
}

constexpr int bases_int[3] = {2, 7, 61};  // intだと、2,7,61で十分
constexpr int bases_ll[7] = {2, 325, 9375, 28178, 450775, 9780504, 1795265022};
template<typename T>
constexpr bool inline is_prime(T n) {
    if (n < 2) {
        return false;
    } else if (n == 2) {
        return true;
    } else if (~n & 1) {
        return false;
    }
    if (n < (1ul << 31)) {
        return miller_rabin<T, dynamic_modint<u32>, bases_int, 3>(n);
    } else {
        return miller_rabin<T, dynamic_modint<u64>, bases_ll, 7>(n);
    }
}

};  // namespace miller
};  // namespace kyopro
#line 3 "test/yuki/No3030.test.cpp"
int main(){
    int n;
    scanf("%d", &n);
    for (int i = 0; i < n; ++i){
        long long x;
        scanf("%lld", &x);
        printf("%lld %c\n", x, kyopro::miller::is_prime(x) ? '1' : '0');
    }
}
yukicoder

結果

ソースコード
