結果

問題 No.3030 ミラー・ラビン素数判定法のテスト
ユーザー nonamaenonamae
提出日時 2022-07-25 08:36:47
言語 C++17
(gcc 12.3.0 + boost 1.83.0)
結果
RE  
実行時間 -
コード長 5,056 bytes
コンパイル時間 1,733 ms
コンパイル使用メモリ 205,128 KB
実行使用メモリ 6,944 KB
最終ジャッジ日時 2024-07-07 05:26:49
合計ジャッジ時間 3,424 ms
ジャッジサーバーID
（参考情報） 		judge5 / judge2
このコードへのチャレンジ
（要ログイン）

テストケース

テストケース表示
入力 結果 実行時間
実行使用メモリ
testcase_00 AC 1 ms
6,812 KB
testcase_01 AC 2 ms
6,944 KB
testcase_02 AC 1 ms
6,944 KB
testcase_03 AC 1 ms
6,944 KB
testcase_04 RE -
testcase_05 RE -
testcase_06 RE -
testcase_07 RE -
testcase_08 RE -
testcase_09 RE -
権限があれば一括ダウンロードができます

ソースコード

diff # 

#include <bits/stdc++.h>

using i8 = std::int8_t;
using i16 = std::int16_t;
using i32 = std::int32_t;
using i64 = std::int64_t;
using i128 = __int128_t;
using u8 = std::uint8_t;
using u16 = std::uint16_t;
using u32 = std::uint32_t;
using u64 = std::uint64_t;
using u128 = __uint128_t;

template<typename T> using vec = std::vector<T>;
template<typename T> using vvec = std::vector<std::vector<T>>;
template<typename T> using vvvec = std::vector<std::vector<std::vector<T>>>;


int jacobi(i64 a, u64 n) {
    u64 t;
    int j = 1;
    while (a) {
        if (a < 0) {
            a = -a;
            if ((n & 3) == 3) j = -j;
        }
        int s = __builtin_ctzll(a);
        a >>= s;
        if (((n & 7) == 3 || (n & 7) == 5) && (s & 1)) j = -j;
        if ((a & n & 3) == 3) j = -j;
        t = a, a = n, n = t;
        a %= n;
        if (((u64)a) > (n / 2)) a -= n;
    }
    return n == 1 ? j : 0;
}

struct m64 {

    // r * m === -1 (mod 1<<64), n2 === 1<<<128 (mod m)
    inline static u64 m, r, n2;

    static void set_mod(u64 m) {
        assert(m < (1ull << 62));
        assert((m & 1) == 1);
        m64::m = m;
        n2 = -u128(m) % m;
        r = m;
        for (int _ = 0; _ < 5; ++_) r *= 2 - m * r;
        r = -r;
        assert(r * m == -1ull);
    }

    static u64 reduce(u128 b) {
        return (b + u128(u64(b) * r) * m) >> 64;
    }

    u64 x;

    m64() : x(0) {}
    m64(u64 x) : x(reduce(u128(x) * n2)) {}
    u64 val() const {
        u64 y = reduce(x);
        return y >= m ? y - m : y;
    }
    m64 &operator+=(m64 y) {
        x += y.x - (m << 1);
        x = (i64(x) < 0 ? x + (m << 1) : x);
        return *this;
    }
    m64 &operator-=(m64 y) {
        x -= y.x;
        x = (i64(x) < 0 ? x + (m << 1) : x);
        return *this;
    }
    m64 &operator*=(m64 y) {
        x = reduce(u128(x) * y.x);
        return *this;
    }
    m64 &operator/=(m64 y) {
        *this *= y.inv();
        return *this;
    }
    m64 &operator<<=(int y) {
        m64 a(1ull << y);
        *this *= a;
        return *this;
    }
    m64 &operator>>=(int y) {
        m64 a(1ull << y);
        *this /= a;
        return *this;
    }
    m64 operator+(m64 y) const { return m64(*this) += y; }
    m64 operator-(m64 y) const { return m64(*this) -= y; }
    m64 operator*(m64 y) const { return m64(*this) *= y; }
    m64 operator/(m64 y) const { return m64(*this) /= y; }
    m64 operator<<(int y) const { return m64(*this) <<= y; }
    m64 operator>>(int y) const { return m64(*this) >>= y; }
    bool operator==(m64 y) const { return (x >= m ? x-m : x) == (y.x >= m ? y.x-m : y.x); }
    bool operator!=(m64 y) const { return not operator==(y); }
    m64 pow(u64 n) const {
        m64 y = 1, z = *this;
        for ( ; n; n >>= 1, z *= z) if (n & 1) y *= z;
        return y;
    }
    m64 inv() const {
        return this->pow(m - 2);
    }
};

bool is_prime(u64 n) {

    {// [1]
        /** 試し割り法 前処理として偶数と小さな数を判定しておく **/
        if (n == 2 || n == 3 || n == 5 || n == 7) return true;
        if (n % 2 == 0 || n % 3 == 0 || n % 5 == 0 || n % 7 == 0) return false;
        if (n < 121) return n > 1;
    }
    m64::set_mod(n);
    {// [2]
        /** 2を底とするMiller-Rabinテストにより、強フェルマー擬素数を見つける **/
        u64 d = (n - 1) << __builtin_clzll(n - 1);
        m64 t(2);
        for (d <<= 1; d; d <<= 1) {
            t *= t;
            if (d >> 63) t <<= 1;
        }
        if (t != m64(1)) {
            u64 x = (n - 1) & -(n - 1);
            m64 rev(n - 1);
            for (x >>= 1; t != rev; x >>= 1) {
                if (x == 0) return false;
                t *= t;
            }
        }
    }
    {// [3]
        /** 初期値・係数の決定 **/
        i64 D = 5;
        for (int i = 0; jacobi(D, n) != -1 && i < 64; ++i) {
            if (i == 32) {
                u32 k = round(sqrtl(n));
                if (k * k == n) return false;
            }
            if (i & 1) D -= 2;
            else       D += 2;
            D = -D;
        }
        m64 Q(D < 0 ? (1 - D) / 4 % n : n - (D - 1) / 4 % n);
        m64 u(1), v(1), Qn(D < 0 ? (1 - D) / 4 % n : n - (D - 1) / 4 % n);
        m64 D_(D < 0 ? (n + D) % n : D % n );
        u64 k = (n + 1) << __builtin_clzll(n + 1);
        for (k <<= 1; k; k <<= 1) {
            u *= v;
            v = v * v - Qn - Qn;
            Qn *= Qn;
            if (k >> 63) {
                m64 uu = (u + v) >> 1;
                v += D_ * u;
                v >>= 1;
                u = uu;
                Qn *= Q;
            }
        }
        if (u == m64(0) || v == m64(0)) return true;
        u64 x = (n + 1) & ~n;
        for (x >>= 1; x; x >>= 1) {
            u *= v;
            v = v * v - Qn - Qn;
            if (v == m64(0)) return true;
            Qn *= Qn;
        }
    }
    return false;
}

int main(void) {
    int n; scanf("%d", &n);
    while (n--) {
        u64 x; scanf("%lu", &x);
        printf("%lu %d\n", x, is_prime(x));
    }
}
0