結果

問題 No.8030 ミラー・ラビン素数判定法のテスト
ユーザー rusa
提出日時 2018-06-01 15:57:51
言語 C++17
(gcc 13.3.0 + boost 1.87.0)
結果
AC  
実行時間 30 ms / 9,973 ms
コード長 6,591 bytes
コンパイル時間 432 ms
コンパイル使用メモリ 55,720 KB
最終ジャッジ日時 2025-01-05 10:48:06
ジャッジサーバーID
（参考情報） 		judge2 / judge3
このコードへのチャレンジ
（要ログイン）
ファイルパターン 結果
other AC * 10
権限があれば一括ダウンロードができます
コンパイルメッセージ
main.cpp: In function ‘int main()’:
main.cpp:232:8: warning: ignoring return value of ‘int scanf(const char*, ...)’ declared with attribute ‘warn_unused_result’ [-Wunused-result]
  232 |   scanf("%d", &T);
      |   ~~~~~^~~~~~~~~~
main.cpp:235:10: warning: ignoring return value of ‘int scanf(const char*, ...)’ declared with attribute ‘warn_unused_result’ [-Wunused-result]
  235 |     scanf("%llu", &n);
      |     ~~~~~^~~~~~~~~~~~

ソースコード

diff # 

#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 - 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 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)));
  }
}
0