結果

問題 No.8030 ミラー・ラビン素数判定法のテスト
ユーザー pyranine
提出日時 2021-01-02 09:18:08
言語 C++14
(gcc 13.3.0 + boost 1.87.0)
結果
AC  
実行時間 48 ms / 9,973 ms
コード長 6,908 bytes
コンパイル時間 1,564 ms
コンパイル使用メモリ 168,448 KB
実行使用メモリ 5,248 KB
最終ジャッジ日時 2024-11-16 23:36:22
合計ジャッジ時間 2,392 ms
ジャッジサーバーID
（参考情報） 		judge2 / judge3
このコードへのチャレンジ
（要ログイン）
ファイルパターン 結果
other AC * 10
権限があれば一括ダウンロードができます

ソースコード

diff # 

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

using Int8 = int8_t;
using Int16 = int16_t;
using Int32 = int32_t;
using Int64 = int64_t;
using Int128 = __int128_t;

using Word8 = uint8_t;
using Word16 = uint16_t;
using Word32 = uint32_t;
using Word64 = uint64_t;
using Word128 = __uint128_t;

using Int = int_fast64_t;
using Word = uint_fast64_t;

using F32 = float;
using F64 = double;
using F80 = long double;

using VS = vector<string>;
using VVS = vector<vector<string>>;
using VB = vector<bool>;
using VVB = vector<vector<bool>>;
using VI = vector<Int>;
using VW = vector<Word>;
using VVI = vector<vector<Int>>;
using VVW = vector<vector<Word>>;
using PII = pair<Int,Int>;
using PWW = pair<Word,Word>;
using VPII = vector<pair<Int,Int>>;
using VPWW = vector<pair<Word,Word>>;

#define LOOP(n) for(Int _ipiewnsjiw=0; _ipiewnsjiw<(n); _ipiewnsjiw++)
#define REP(i,n) for(Int i=0, i##_len=(n); i<i##_len; ++i)
#define RANGE(i,a,b) for(Int i=(a), i##_len(b); i<=i##_len; ++i)
#define SZ(obj) ((Int)(obj).size())
#define UNIQUE(obj) (obj).erase(unique((obj).begin(),(obj).end()),(obj).end())
#define ALL(obj) (obj).begin(),(obj).end()
#define RALL(obj) (obj).rbegin(),(obj).rend()

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(Word64 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<Word64, Word128, Int64>;
using UnsafeMod32 = UnsafeMod<Word32, Word64, int>;
template <> Word64 UnsafeMod64::mod = 0;
template <> Word64 UnsafeMod64::inv = 0;
template <> Word64 UnsafeMod64::r2 = 0;
template <> Word32 UnsafeMod32::mod = 0;
template <> Word32 UnsafeMod32::inv = 0;
template <> Word32 UnsafeMod32::r2 = 0;

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(Word64 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<Word64, Word128, Int64>;
using Mod32 = Mod<Word32, Word64, int>;
template <> Word64 Mod64::mod = 0;
template <> Word64 Mod64::inv = 0;
template <> Word64 Mod64::r2 = 0;
template <> Word64 Mod64::mask = 0;
template <> Word64 Mod64::one = 0;
template <> int Mod64::shift = 0;
template <> Word32 Mod32::mod = 0;
template <> Word32 Mod32::inv = 0;
template <> Word32 Mod32::r2 = 0;
template <> Word32 Mod32::mask = 0;
template <> Word32 Mod32::one = 0;
template <> int Mod32::shift = 0;


template <class word, class mod>
bool is_prime(word n, const Word32* bases, int m) {
  mod::set_mod(n);
  int s = __builtin_ctzll(n - 1);
  word d = (n - 1) >> s;
  mod one{1}, minus_one{n - 1};
  for (int i = 0, j; i < m; ++i) {
    mod a = mod(bases[i]).pow(d);
    if (a == one || a == minus_one) continue;
    for (j = s - 1; j > 0; --j) {
      if ((a *= a) == minus_one) break;
    }
    if (j == 0) return false;
  }
  return true;
}

bool miller_rabin(Word n) {
  static const Word32 bases[][7] = {
    {2, 3},
    {2, 299417},
    {2, 7, 61},
    {15, 176006322, Word32(4221622697)},
    {2, 2570940, 211991001, Word32(3749873356)},
    {2, 2570940, 880937, 610386380, Word32(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 < (Word32(1) << 31)) return is_prime<Word32, UnsafeMod32>(n, bases[x], y);
  else if (n < (Word64(1) << 63)) return is_prime<Word64, UnsafeMod64>(n, bases[x], y);
  else assert(0);
}


int main() {
  int Q; cin >> Q;
  LOOP(Q) {
    Word x; cin >> x;
    cout << x << ' ' << miller_rabin(x) << '\n';
  }
  return 0;
}
0