結果
問題 | No.8030 ミラー・ラビン素数判定法のテスト |
ユーザー | zimpha |
提出日時 | 2018-04-21 15:36:27 |
言語 | C++17 (gcc 12.3.0 + boost 1.83.0) |
結果 |
WA
|
実行時間 | - |
コード長 | 6,586 bytes |
コンパイル時間 | 624 ms |
コンパイル使用メモリ | 53,504 KB |
実行使用メモリ | 6,824 KB |
最終ジャッジ日時 | 2024-11-18 16:26:59 |
合計ジャッジ時間 | 1,628 ms |
ジャッジサーバーID (参考情報) |
judge2 / judge5 |
(要ログイン)
テストケース
テストケース表示入力 | 結果 | 実行時間 実行使用メモリ |
---|---|---|
testcase_00 | WA | - |
testcase_01 | WA | - |
testcase_02 | WA | - |
testcase_03 | WA | - |
testcase_04 | WA | - |
testcase_05 | WA | - |
testcase_06 | WA | - |
testcase_07 | WA | - |
testcase_08 | WA | - |
testcase_09 | WA | - |
ソースコード
#include <cassert> using uint32 = unsigned int; using int64 = long long; using uint64 = unsigned long long; using uint128 = __uint128_t; // Montgomery modular multiplication -- about 7x faster // ensure mod is an odd number, use after call `set_mod` method template<typename word, typename dword, typename sword> struct UnsafeMod { UnsafeMod(): x(0) {} UnsafeMod(word _x): x(init(_x)) {} bool operator == (const UnsafeMod& rhs) const { return x == rhs.x; } bool operator != (const UnsafeMod& rhs) const { return x != rhs.x; } UnsafeMod& operator += (const UnsafeMod& rhs) { if ((x += rhs.x) >= mod) x -= mod; return *this; } UnsafeMod& operator -= (const UnsafeMod& rhs) { if (sword(x -= rhs.x) < 0) x += mod; return *this; } UnsafeMod& operator *= (const UnsafeMod& rhs) { x = reduce(dword(x) * rhs.x); return *this; } UnsafeMod operator + (const UnsafeMod &rhs) const { return UnsafeMod(*this) += rhs; } UnsafeMod operator - (const UnsafeMod &rhs) const { return UnsafeMod(*this) -= rhs; } UnsafeMod operator * (const UnsafeMod &rhs) const { return UnsafeMod(*this) *= rhs; } UnsafeMod pow(uint64 e) const { UnsafeMod ret(1); for (UnsafeMod base = *this; e; e >>= 1, base *= base) { if (e & 1) ret *= base; } return ret; } word get() const { return reduce(x); } static constexpr int word_bits = sizeof(word) * 8; static word modulus() { return mod; } static word init(word w) { return reduce(dword(w) * r2); } static void set_mod(word m) { mod = m; inv = mul_inv(mod); r2 = -dword(mod) % mod; } static word reduce(dword x) { word y = word(x >> word_bits) - word((dword(word(x) * inv) * mod) >> word_bits); return sword(y) < 0 ? y + mod : y; } static word mul_inv(word n, int e = 6, word x = 1) { return !e ? x : mul_inv(n, e - 1, x * (2 - x * n)); } static word mod, inv, r2; word x; }; using UnsafeMod64 = UnsafeMod<uint64, uint128, int64>; using UnsafeMod32 = UnsafeMod<uint32, uint64, int>; template <> uint64 UnsafeMod64::mod = 0; template <> uint64 UnsafeMod64::inv = 0; template <> uint64 UnsafeMod64::r2 = 0; template <> uint32 UnsafeMod32::mod = 0; template <> uint32 UnsafeMod32::inv = 0; template <> uint32 UnsafeMod32::r2 = 0; // in this version mod can be any positive number // speed is about 5x than usual mod template<typename word, typename dword, typename sword> struct Mod { Mod() : x(0) {} Mod(word _x) : x(init(_x)) {} Mod& operator += (const Mod& rhs) { word hi = (x >> shift) + (rhs.x >> shift) - mod; if (sword(hi) < 0) hi += mod; x = hi << shift | ((x + rhs.x) & mask); return *this; } Mod& operator -= (const Mod& rhs) { word hi = (x >> shift) - (rhs.x >> shift); if (sword(hi) < 0) hi += mod; x = hi << shift | ((x - rhs.x) & mask); return *this; } Mod& operator *= (const Mod& rhs) { x = reduce(x, rhs.x); return *this; } Mod operator + (const Mod& rhs) const { return Mod(*this) += rhs; } Mod operator - (const Mod& rhs) const { return Mod(*this) -= rhs; } Mod operator * (const Mod& rhs) const { return Mod(*this) *= rhs; } word get() const { word ret = reduce(x, one); word r1 = ret >> shift; return mod * (((ret - r1) * inv) & mask) + r1; } Mod pow(uint64 e) const { Mod ret = Mod(1); for (Mod base = *this; e; e >>= 1, base *= base) { if (e & 1) ret *= base; } return ret; } static constexpr int word_bits = sizeof(word) * 8; static void set_mod(word m) { shift = __builtin_ctzll(m); mask = (word(1) << shift) - 1; mod = m >> shift; inv = mul_inv(mod); assert(mod * inv == 1); r2 = -dword(mod) % mod; one = word(1) << shift | 1; } static word modulus() { return mod << shift; } static word init(word x) { return reduce_odd(dword(x) * r2) << shift | (x & mask); } static word reduce_odd(dword x) { word y = word(x >> word_bits) - word((dword(word(x) * inv) * mod) >> word_bits); return sword(y) < 0 ? y + mod : y; } static word reduce(word x0, word x1) { word y = reduce_odd(dword(x0 >> shift) * (x1 >> shift)); return y << shift | ((x0 * x1) & mask); } static word mul_inv(word n, int e = 6, word x = 1) { return !e ? x : mul_inv(n, e - 1, x * (2 - x * n)); } static word mod, inv, r2, mask, one; static int shift; word x; }; using Mod64 = Mod<uint64, uint128, int64>; using Mod32 = Mod<uint32, uint64, int>; template <> uint64 Mod64::mod = 0; template <> uint64 Mod64::inv = 0; template <> uint64 Mod64::r2 = 0; template <> uint64 Mod64::mask = 0; template <> uint64 Mod64::one = 0; template <> int Mod64::shift = 0; template <> uint32 Mod32::mod = 0; template <> uint32 Mod32::inv = 0; template <> uint32 Mod32::r2 = 0; template <> uint32 Mod32::mask = 0; template <> uint32 Mod32::one = 0; template <> int Mod32::shift = 0; #include <cstdio> #include <vector> #include <cstdlib> #include <algorithm> namespace primes { template <class word, class mod> bool composite(word n, const uint32* bases, int m) { mod::set_mod(n); int s = __builtin_ctzll(n); word d = (n - 1) >> s; mod one{1}, minus_one{n - 1}; for (int i = 0, j = s - 1; i < m; ++i) { mod a = mod(bases[i]).pow(d); if (a == one || a == minus_one) continue; for (; j > 0; --j) { if ((a *= a) == minus_one) break; } if (j == 0) return true; } return false; } bool is_prime(uint64 n) { // reference: http://miller-rabin.appspot.com static const uint32 bases[][7] = { {2, 3}, {2, 299417}, {2, 7, 61}, {15, 176006322, uint32(4221622697)}, {2, 2570940, 211991001, uint32(3749873356)}, {2, 2570940, 880937, 610386380, uint32(4130785767)}, {2, 325, 9375, 28178, 450775, 9780504, 1795265022} }; if (n <= 1) return false; if (!(n & 1)) return n == 2; if (n <= 8) return true; int x = 6, y = 7; if (n < 1373653) x = 0, y = 2; else if (n < 19471033) x = 1, y = 2; else if (n < 4759123141) x = 2, y = 3; else if (n < 154639673381) x = y = 3; else if (n < 47636622961201) x = y = 4; else if (n < 3770579582154547) x = y = 5; if (n < (uint32(1) << 31)) { return !composite<uint32, UnsafeMod32>(n, bases[x], y); } else if (n < (uint64(1) << 63)) { return !composite<uint64, UnsafeMod64>(n, bases[x], y); } else assert(0); } } int main() { int T; scanf("%d", &T); for (int cas = 1; cas <= T; ++cas) { uint64 n; scanf("%llu", &n); printf("%llu %d\n", n, int(primes::is_prime(n))); } }