結果

問題 No.8030 ミラー・ラビン素数判定法のテスト
ユーザー こまる
提出日時 2020-11-07 15:04:08
言語 C++17
(gcc 13.3.0 + boost 1.87.0)
結果
WA  
実行時間 -
コード長 8,362 bytes
コンパイル時間 4,245 ms
コンパイル使用メモリ 218,460 KB
最終ジャッジ日時 2025-01-15 21:30:29
ジャッジサーバーID
（参考情報） 		judge3 / judge1
このコードへのチャレンジ
（要ログイン）
ファイルパターン 結果
other WA * 10
権限があれば一括ダウンロードができます
コンパイルメッセージ
main.cpp:354:39: warning: multi-character character constant [-Wmultichar]
  354 |     wt(primes::miller_rabin<i64>(x) ? ' 1\n' : ' 0\n');
      |                                       ^~~~~~
main.cpp:354:48: warning: multi-character character constant [-Wmultichar]
  354 |     wt(primes::miller_rabin<i64>(x) ? ' 1\n' : ' 0\n');
      |                                                ^~~~~~

ソースコード

diff # 

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

using u32  = uint_fast32_t;
using u64  = uint_fast64_t;
using u128 = __uint128_t;
using i32  = int_fast32_t;
using i64  = int_fast64_t;
using i128 = __int128_t;

int isqrt(i64 n) {
  return sqrtl(n);
}

namespace primes {
using ll = long long;
struct mod64 {
  u64 n;
  static u64 mod, inv, r2;
  mod64() : n(0) {}
  mod64(u64 x) : n(init(x)) {}
  static u64 init(u64 w) { return reduce(u128(w) * r2); }
  static void set_mod(u64 m) {
    mod = inv = m;
    for (int i = 0; i < 5; ++i) inv *= 2 - inv * m;
    r2 = -u128(m) % m;
  }
  static u64 reduce(u128 x) {
    u64 y = u64(x >> 64) - u64((u128(u64(x) * inv) * mod) >> 64);
    return ll(y) < 0 ? y + mod : y;
  };
  mod64& operator+=(mod64 x) { n += x.n - mod; if(ll(n) < 0) n += mod; return *this; }
  mod64 operator+(mod64 x) const { return mod64(*this) += x; }
  mod64& operator*=(mod64 x) { n = reduce(u128(n) * x.n);  return *this; }
  mod64 operator*(mod64 x) const { return mod64(*this) *= x; }
  u64 val() const { return reduce(n); }
};

u64 mod64::mod, mod64::inv, mod64::r2;

bool suspect(u64 a, u64 s, u64 d, u64 n){
  mod64::set_mod(n);
  mod64 x(1), xx(a), one(x), minusone(n-1);
  while(d > 0){
    if(d&1) x = x * xx;
    xx = xx * xx;
    d >>= 1;
  }
  if (x.n == one.n) return true;
  for (int r = 0; r < s; ++r) {
    if(x.n == minusone.n) return true;
    x = x * x;
  }
  return false;
}

template<class T>
bool miller_rabin(T n){
  if (n <= 1 || (n > 2 && n % 2 == 0)) return false;
  u64 d = n - 1, s = 0;
  while (!(d&1)) {++s; d >>= 1;}
  static const u64 v[] = {2, 325, 9375, 28178, 450775, 9780504, 1795265022};
  static const u64 v_small[] = {2, 7, 61};
  if(n <= 4759123141LL){
    for (auto &&p : v_small) {
      if(p >= n) break;
      if(!suspect(p, s, d, n)) return false;
    }
  }else {
    for (auto &&p : v) {
      if(p >= n) break;
      if(!suspect(p, s, d, n)) return false;
    }
  }
  return true;
}

template<typename T>
struct ExactDiv {
  T t, i, val;
  ExactDiv() {}
  ExactDiv(T n) : t(T(-1) / n), i(mul_inv(n)) , val(n) {};
  T mul_inv(T n) {
    T x = n;
    for (int i = 0; i < 5; ++i) x *= 2 - n * x;
    return x;
  }
  bool divide(T n) const {
    if(val == 2) return !(n & 1);
    return n * this->i <= this->t;
  }
};

vector<ExactDiv<u64>> get_prime(int n){
  if(n <= 1) return vector<ExactDiv<u64>>();
  vector<bool> is_prime(n+1, true);
  vector<ExactDiv<u64>> prime;
  is_prime[0] = is_prime[1] = false;
  for (int i = 2; i <= n; ++i) {
    if(is_prime[i]) prime.emplace_back(i);
    for (auto &&j : prime){
      if(i*j.val > n) break;
      is_prime[i*j.val] = false;
      if(j.divide(i)) break;
    }
  }
  return prime;
}
const auto primes = get_prime(50000);

random_device rng;


template<class T>
T pollard_rho2(T n) {
  uniform_int_distribution<T> ra(1, n-1);
  mod64::set_mod(n);
  while(true){
    u64 c_ = ra(rng), g = 1, r = 1, m = 500;
    while(c_ == n-2) c_ = ra(rng);
    mod64 y(ra(rng)), xx(0), c(c_), ys(0), q(1);
    while(g == 1){
      xx.n = y.n;
      for (int i = 1; i <= r; ++i) {
        y *= y; y += c;
      }
      T k = 0; g = 1;
      while(k < r && g == 1){
        for (int i = 1; i <= (m > (r-k) ? (r-k) : m); ++i) {
          ys.n = y.n;
          y *= y; y += c;
          u64 xxx = xx.val(), yyy = y.val();
          q *= mod64(xxx > yyy ? xxx - yyy : yyy - xxx);
        }
        g = __gcd<ll>(q.val(), n);
        k += m;
      }
      r *= 2;
    }
    if(g == n) g = 1;
    while (g == 1){
      ys *= ys; ys += c;
      u64 xxx = xx.val(), yyy = ys.val();
      g = __gcd<ll>(xxx > yyy ? xxx - yyy : yyy - xxx, n);
    }
    if (g != n && miller_rabin(g)) return g;
  }
}

template<class T>
vector<T> prime_factor(T n, int d = 0){
  vector<T> a, res;
  if(!d) for (auto &&i : primes) {
      while (i.divide(n)){
        res.emplace_back(i.val);
        n /= i.val;
      }
    }
  while(n != 1){
    if(miller_rabin(n)){
      a.emplace_back(n);
      break;
    }
    T x = pollard_rho2(n);
    n /= x;
    a.emplace_back(x);
  }
  for (auto &&i : a) {
    if (miller_rabin(i)) {
      res.emplace_back(i);
    } else {
      vector<T> b = prime_factor(i, d + 1);
      for (auto &&j : b) res.emplace_back(j);
    }
  }
  sort(res.begin(),res.end());
  return res;
}

__attribute__((target("avx"), optimize("O3", "unroll-loops")))
i64 prime_pi(const i64 N) {
  if (N <= 1) return 0;
  if (N == 2) return 1;
  const int v = isqrt(N);
  int s = (v + 1) / 2;
  vector<int> smalls(s); for (int i = 1; i < s; ++i) smalls[i] = i;
  vector<int> roughs(s); for (int i = 0; i < s; ++i) roughs[i] = 2 * i + 1;
  vector<i64> larges(s); for (int i = 0; i < s; ++i) larges[i] = (N / (2 * i + 1) - 1) / 2;
  vector<bool> skip(v + 1);
  const auto divide = [] (i64 n, i64 d) -> int { return double(n) / d; };
  const auto half = [] (int n) -> int { return (n - 1) >> 1; };
  int pc = 0;
  for (int p = 3; p <= v; p += 2) if (!skip[p]) {
    int q = p * p;
    if (i64(q) * q > N) break;
    skip[p] = true;
    for (int i = q; i <= v; i += 2 * p) skip[i] = true;
    int ns = 0;
    for (int k = 0; k < s; ++k) {
      int i = roughs[k];
      if (skip[i]) continue;
      i64 d = i64(i) * p;
      larges[ns] = larges[k] - (d <= v ? larges[smalls[d >> 1] - pc] : smalls[half(divide(N, d))]) + pc;
      roughs[ns++] = i;
    }
    s = ns;
    for (int i = half(v), j = ((v / p) - 1) | 1; j >= p; j -= 2) {
      int c = smalls[j >> 1] - pc;
      for (int e = (j * p) >> 1; i >= e; --i) smalls[i] -= c;
    }
    ++pc;
  }
  larges[0] += i64(s + 2 * (pc - 1)) * (s - 1) / 2;
  for (int k = 1; k < s; ++k) larges[0] -= larges[k];
  for (int l = 1; l < s; ++l) {
    int q = roughs[l];
    i64 M = N / q;
    int e = smalls[half(M / q)] - pc;
    if (e < l + 1) break;
    i64 t = 0;
    for (int k = l + 1; k <= e; ++k) t += smalls[half(divide(M, roughs[k]))];
    larges[0] += t - i64(e - l) * (pc + l - 1);
  }
  return larges[0] + 1;
}


} // namespace primes

namespace fastio {
static constexpr int SZ = 1 << 17;
char ibuf[SZ], obuf[SZ];
int pil = 0, pir = 0, por = 0;

struct Pre {
  char num[40000];
  constexpr Pre() : num() {
    for (int i = 0; i < 10000; i++) {
      int n = i;
      for (int j = 3; j >= 0; j--) {
        num[i * 4 + j] = n % 10 + '0';
        n /= 10;
      }
    }
  }
} constexpr pre;

inline void load() {
  memcpy(ibuf, ibuf + pil, pir - pil);
  pir = pir - pil + fread(ibuf + pir - pil, 1, SZ - pir + pil, stdin);
  pil = 0;
}
inline void flush() {
  fwrite(obuf, 1, por, stdout);
  por = 0;
}

inline void rd(char& c) { c = ibuf[pil++]; }
template <typename T>
inline void rd(T& x) {
  if (pil + 32 > pir) load();
  char c;
  do
    c = ibuf[pil++];
  while (c < '-');
  bool minus = 0;
  if (c == '-') {
    minus = 1;
    c = ibuf[pil++];
  }
  x = 0;
  while (c >= '0') {
    x = x * 10 + (c & 15);
    c = ibuf[pil++];
  }
  if (minus) x = -x;
}
inline void rd() {}
template <typename Head, typename... Tail>
inline void rd(Head& head, Tail&... tail) {
  rd(head);
  rd(tail...);
}

inline void wt(char c) { obuf[por++] = c; }
template <typename T>
inline void wt(T x) {
  if (por > SZ - 32) flush();
  if (!x) {
    obuf[por++] = '0';
    return;
  }
  if (x < 0) {
    obuf[por++] = '-';
    x = -x;
  }
  int i = 12;
  char buf[16];
  while (x >= 10000) {
    memcpy(buf + i, pre.num + (x % 10000) * 4, 4);
    x /= 10000;
    i -= 4;
  }
  if (x < 100) {
    if (x < 10) {
      wt(char('0' + char(x)));
    } else {
      uint32_t q = (uint32_t(x) * 205) >> 11;
      uint32_t r = uint32_t(x) - q * 10;
      obuf[por + 0] = '0' + q;
      obuf[por + 1] = '0' + r;
      por += 2;
    }
  } else {
    if (x < 1000) {
      memcpy(obuf + por, pre.num + (x << 2) + 1, 3);
      por += 3;
    } else {
      memcpy(obuf + por, pre.num + (x << 2), 4);
      por += 4;
    }
  }
  memcpy(obuf + por, buf + i + 4, 12 - i);
  por += 12 - i;
}

inline void wt() {}
template <typename Head, typename... Tail>
inline void wt(Head head, Tail... tail) {
  wt(head);
  wt(tail...);
}
template <typename T>
inline void wtn(T x) {
  wt(x, '\n');
}

struct Dummy {
  Dummy() { atexit(flush); }
} dummy;

}  // namespace fastio
using fastio::rd;
using fastio::wt;
using fastio::wtn;

int main() {
  int n;
  rd(n);
  while(n--) {
    i64 x;
    rd(x);
    wt(x);
    wt(primes::miller_rabin<i64>(x) ? ' 1\n' : ' 0\n');
  }
  return 0;
}
0