#560010 (C++17) No.3030 ミラー・ラビン素数判定法のテスト

提出ソース

結果

問題	No.3030 ミラー・ラビン素数判定法のテスト
ユーザー	Pyこまる
提出日時	2020-09-29 08:10:57
言語	C++17 (gcc 12.3.0 + boost 1.83.0)
結果	RE
実行時間	-
コード長	7,278 bytes
コンパイル時間	859 ms
コンパイル使用メモリ	77,592 KB
実行使用メモリ	6,824 KB
最終ジャッジ日時	2024-11-18 18:27:59
合計ジャッジ時間	2,121 ms
ジャッジサーバーID （参考情報）	judge5 / judge4

このコードへのチャレンジ
（要ログイン）

テストケース

テストケース表示

入力	結果	実行時間実行使用メモリ
testcase_00	AC	2 ms 6,816 KB
testcase_01	AC	2 ms 6,816 KB
testcase_02	AC	2 ms 6,816 KB
testcase_03	AC	2 ms 6,820 KB
testcase_04	RE	-
testcase_05	RE	-
testcase_06	RE	-
testcase_07	RE	-
testcase_08	RE	-
testcase_09	RE	-

権限があれば一括ダウンロードができます

ソースコード

raw source code

#include <iostream>
#include <vector>
#include <cassert>
using namespace std;

using i32 = int32_t;
using i64 = int64_t;
using i128 = __int128_t;
using u32 = uint32_t;
using u64 = uint64_t;
using u128 = __uint128_t;

template <typename T = i32>
T gcd(T a, T b) {
    while (b) swap(a %= b, b);
    return a;
}

u64 gcd_fast(u64 x, u64 y) {
    if (x == 0 || y == 0) return x + y;
    u32 bx = __builtin_ctz(x);
    u32 by = __builtin_ctz(y);
    u32 k = min(bx, by);
    x >>= bx, y >>= by;
    while (x != 0 && y != 0)
        if (x == y) return x << k;
        else if (x > y) x = (x - y) >> __builtin_ctz(x - y);
        else y = (y - x) >> __builtin_ctz(y - x);
    return (x + y) << k;
}

template <typename T = i32>
T inv(T a, T p) {
    T b = p, x = 1, y = 0;
    while (a) {
        T q = b % a;
        swap(a, b /= a);
        swap(x, y -= q * x);
    }
    assert(b == 1);
    return y < 0 ? y + p : y;
}

template <typename T = int32_t, typename U = int64_t>
T modpow(T a, U n, T p) {
    T ret = 1;
    for (; n; n >>= 1, a = U(a) * a % p)
        if (n & 1) ret = U(ret) * a % p;
    return ret;
}

struct ArbitraryLazyMontgomeryModInt {
    using mint = ArbitraryLazyMontgomeryModInt;
    using i32 = int32_t;
    using u32 = uint32_t;
    using u64 = uint64_t;

    static u32 mod;
    static u32 r;
    static u32 n2;

    static u32 get_r() {
        u32 ret = mod;
        for (i32 i = 0; i < 4; ++i) ret *= 2 - mod * ret;
        return ret;
    }

    static void set_mod(u32 m) {
        assert(m < (1 << 30));
        assert((m & 1) == 1);
        mod = m;
        n2 = -u64(m) % m;
        r = get_r();
        assert(r * mod == 1);
    }

    u32 a;

    ArbitraryLazyMontgomeryModInt() : a(0) {}
    ArbitraryLazyMontgomeryModInt(const int64_t &b) : a(reduce(u64(b % mod + mod) * n2)){};

    static u32 reduce(const u64 &b) {
        return (b + u64(u32(b) * u32(-r)) * mod) >> 32;
    }

    mint &operator+=(const mint &b) {
        if (i32(a += b.a - 2 * mod) < 0) a += 2 * mod;
        return *this;
    }

    mint &operator-=(const mint &b) {
        if (i32(a -= b.a) < 0) a += 2 * mod;
        return *this;
    }

    mint &operator*=(const mint &b) {
        a = reduce(u64(a) * b.a);
        return *this;
    }

    mint &operator/=(const mint &b) {
        *this *= b.inverse();
        return *this;
    }

    mint operator+(const mint &b) const { return mint(*this) += b; }
    mint operator-(const mint &b) const { return mint(*this) -= b; }
    mint operator*(const mint &b) const { return mint(*this) *= b; }
    mint operator/(const mint &b) const { return mint(*this) /= b; }
    bool operator==(const mint &b) const {
        return (a >= mod ? a - mod : a) == (b.a >= mod ? b.a - mod : b.a);
    }
    bool operator!=(const mint &b) const {
        return (a >= mod ? a - mod : a) != (b.a >= mod ? b.a - mod : b.a);
    }
    mint operator-() const { return mint() - mint(*this); }

    mint pow(u64 n) const {
        mint ret(1), mul(*this);
        while (n > 0) {
            if (n & 1) ret *= mul;
            mul *= mul;
            n >>= 1;
        }
        return ret;
    }

    friend ostream &operator<<(ostream &os, const mint &b) {
        return os << b.get();
    }

    friend istream &operator>>(istream &is, mint &b) {
        int64_t t;
        is >> t;
        b = ArbitraryLazyMontgomeryModInt(t);
        return (is);
    }

    mint inverse() const { return pow(mod - 2); }

    u32 get() const {
        u32 ret = reduce(a);
        return ret >= mod ? ret - mod : ret;
    }

    static u32 get_mod() { return mod; }
};
typename ArbitraryLazyMontgomeryModInt::u32 ArbitraryLazyMontgomeryModInt::mod;
typename ArbitraryLazyMontgomeryModInt::u32 ArbitraryLazyMontgomeryModInt::r;
typename ArbitraryLazyMontgomeryModInt::u32 ArbitraryLazyMontgomeryModInt::n2;


struct montgomery64 {
    using mint = montgomery64;
    using i64 = int64_t;
    using u64 = uint64_t;
    using u128 = __uint128_t;

    static u64 mod;
    static u64 r;
    static u64 n2;

    static u64 get_r() {
        u64 ret = mod;
        for (i64 i = 0; i < 5; ++i) ret *= 2 - mod * ret;
        return ret;
    }

    static void set_mod(u64 m) {
        assert(m < (1LL << 62));
        assert((m & 1) == 1);
        mod = m;
        n2 = -u128(m) % m;
        r = get_r();
        assert(r * mod == 1);
    }

    u64 a;

    montgomery64() : a(0) {}
    montgomery64(const int64_t &b) : a(reduce((u128(b) + mod) * n2)){};

    static u64 reduce(const u128 &b) {
      return (b + u128(u64(b) * u64(-r)) * mod) >> 64;
    }

    mint &operator+=(const mint &b) {
        if (i64(a += b.a - 2 * mod) < 0) a += 2 * mod;
        return *this;
    }

    mint &operator-=(const mint &b) {
        if (i64(a -= b.a) < 0) a += 2 * mod;
        return *this;
    }

    mint &operator*=(const mint &b) {
        a = reduce(u128(a) * b.a);
        return *this;
    }

    mint &operator/=(const mint &b) {
        *this *= b.inverse();
        return *this;
    }

    mint operator+(const mint &b) const { return mint(*this) += b; }
    mint operator-(const mint &b) const { return mint(*this) -= b; }
    mint operator*(const mint &b) const { return mint(*this) *= b; }
    mint operator/(const mint &b) const { return mint(*this) /= b; }
    bool operator==(const mint &b) const {
        return (a >= mod ? a - mod : a) == (b.a >= mod ? b.a - mod : b.a);
    }
    bool operator!=(const mint &b) const {
        return (a >= mod ? a - mod : a) != (b.a >= mod ? b.a - mod : b.a);
    }
    mint operator-() const { return mint() - mint(*this); }

    mint pow(u128 n) const {
        mint ret(1), mul(*this);
        while (n > 0) {
            if (n & 1) ret *= mul;
            mul *= mul;
            n >>= 1;
        }
        return ret;
    }

    friend ostream &operator<<(ostream &os, const mint &b) {
        return os << b.get();
    }

    friend istream &operator>>(istream &is, mint &b) {
        int64_t t;
        is >> t;
        b = montgomery64(t);
        return (is);
    }

      mint inverse() const { return pow(mod - 2); }

    u64 get() const {
        u64 ret = reduce(a);
        return ret >= mod ? ret - mod : ret;
    }

    static u64 get_mod() { return mod; }
};
typename montgomery64::u64 montgomery64::mod, montgomery64::r, montgomery64::n2;

template <typename mint>
bool miller_rabin(u64 n, vector<u64> as) {
    if (mint::get_mod() != n) mint::set_mod(n);
    u64 d = n - 1;
    while (~d&1) d >>= 1;
    mint e{1}, rev{i64(n-1)};
    for (u64 a : as) {
        if (n <= a) break;
        u64 t = d;
        mint y = mint(a).pow(t);
        while(t != n - 1 && y != e && y != rev) {
            y *= y;
            t *= 2;
        }
        if (y != rev && t % 2 == 0) return false;
    }
    return true;
}

bool is_prime(u64 n) {
    if (n <= 3) return n == 2 || n == 3;
    if (n % 2 == 0) return false;
    if (n < (1LL << 30)) return miller_rabin<ArbitraryLazyMontgomeryModInt>(n, {2, 7, 61});
    else return miller_rabin<montgomery64>(n, {2, 325, 9375, 28178, 450775, 9780504, 1795265022});
}


int main() {
    int query;
    u64 x;
    cin >> query;
    for (int i = 0; i < query; i++) {
        cin >> x;
        cout << x << " " << (is_prime(x) ? 1 : 0) << endl;
    }
    return 0;
}

yukicoder

結果

テストケース

ソースコード