結果

問題 No.8030 ミラー・ラビン素数判定法のテスト
ユーザー satonakatakumi
提出日時 2021-08-23 16:53:18
言語 C
(gcc 13.3.0)
結果
WA  
実行時間 -
コード長 3,298 bytes
コンパイル時間 1,339 ms
コンパイル使用メモリ 33,152 KB
実行使用メモリ 5,248 KB
最終ジャッジ日時 2024-11-07 19:32:59
合計ジャッジ時間 1,826 ms
ジャッジサーバーID
（参考情報） 		judge1 / judge4
このコードへのチャレンジ
（要ログイン）
ファイルパターン 結果
other AC * 4 WA * 6
権限があれば一括ダウンロードができます

ソースコード

diff # 

#include <assert.h>
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <stdbool.h>
#include <time.h>
#include <math.h>

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

static inline u64 ctz(u64 n) { return __builtin_ctzll(n); }
static inline u64 clz(u64 n) { return __builtin_clzll(n); }
static inline u64 popcnt(u64 n) { return __builtin_popcountll(n); }
static inline u64 gcd(u64 u, u64 v) {
  if (u == 0) return v;
  if (v == 0) return u;
  u64 k;
  for (k = 0; ((u | v) & 1) == 0; ++k) {
    u >>= 1;
    v >>= 1;
  }
  while ((u & 1) == 0) u >>= 1;
  do {
    while ((v & 1) == 0) v >>= 1;
    if (u > v) {
      u ^= v;
      v ^= u;
      u ^= v;
    }
    v = v - u;
  } while (v != 0);
  return u << k;
}

typedef struct Montgomery_t {
  u64 m, r, n2;
  u64 x;
} m64;

m64 make_montgomery(u64 mod) {
  m64 result;
  assert(mod != 0);
  assert(mod & 1 != 0);
  result.m = mod;
  result.n2 = -(u128)(mod) % mod;
  result.r = mod;
  for (int _ = 0; _ < 5; _++)
    result.r *= 2 - result.m * result.r;
  result.r = -result.r;
  return result;
}
u64 MR(m64 p, u128 b) { return (b + (u128)((u64)b * p.r) * p.m) >> 64; }
u64 RM(m64 p) {
  u64 y = MR(p, p.x);
  return y >= p.m ? y - p.m : y;
}
m64 construct(u64 z, u64 mod) {
  m64 p = make_montgomery(mod);
  p.x = MR(p, (u128)z * p.n2);
  return p;
}
m64 plus(m64 p, m64 q) {
  p.x += q.x - (q.m << 1);
  p.x = ((i64)p.x < 0 ? p.x + (p.m << 1) : p.x);
  return p;
}
m64 minus(m64 p, m64 q) {
  p.x -= q.x;
  p.x = ((i64)p.x < 0 ? p.x + (p.m << 1) : p.x);
  return p;
}
m64 mult(m64 p, m64 q) {
  p.x = MR(p, (u128)p.x * q.x);
  return p;
}
m64 powmod(m64 p, u64 n) {
  m64 y = construct(1, p.m);
  m64 z = p;
  for (; n; n >>= 1, z = mult(z, z))
    if (n & 1)
      y = mult(y, z);
  return y;
}
bool eq(m64 p, m64 q) { return (p.x >= p.m ? p.x - p.m : p.x) == (q.x >= q.m ? q.x - q.m : q.x); }
bool not_eq(m64 p, m64 q) { return !eq(p, q); }

bool composite(u64 n, const u32* bases, int m) {
  u64 s = ctz(n - 1);
  u64 d = (n - 1) >> s;
  for (int i = 0, j; i < m; ++i) {
    m64 a = powmod(construct(bases[i], n), d);
    if (1 == RM(a) || RM(a) == n - 1) continue;
    for (j = s - 1; j > 0; --j) {
      a = mult(a, a);
      if (RM(a) == n - 1) break;
    }
    if (j == 0) return true;
  }
  return false;
}

bool is_prime(u64 n) {
  static const u32 bases[][7] = {
    {2, 3},
    {2, 299417},
    {2, 7, 61},
    {15, 176006322, (u32)(4221622697)},
    {2, 2570940, 211991001, (u32)(3749873356)},
    {2, 2570940, 880937, 610386380, (u32)(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;
  return !composite(n, bases[x], y);
}

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