結果

問題 No.3030 ミラー・ラビン素数判定法のテスト
ユーザー nonamaenonamae
提出日時 2021-11-07 00:23:58
言語 C
(gcc 12.3.0)
結果
AC  
実行時間 16 ms / 9,973 ms
コード長 3,866 bytes
コンパイル時間 257 ms
コンパイル使用メモリ 33,552 KB
実行使用メモリ 4,380 KB
最終ジャッジ日時 2023-08-10 16:33:56
合計ジャッジ時間 978 ms
ジャッジサーバーID
（参考情報） 		judge12 / judge15
このコードへのチャレンジ（β）

テストケース

テストケース表示
入力 結果 実行時間
実行使用メモリ
testcase_00 AC 1 ms
4,376 KB
testcase_01 AC 0 ms
4,376 KB
testcase_02 AC 0 ms
4,376 KB
testcase_03 AC 1 ms
4,380 KB
testcase_04 AC 13 ms
4,376 KB
testcase_05 AC 14 ms
4,376 KB
testcase_06 AC 11 ms
4,380 KB
testcase_07 AC 11 ms
4,380 KB
testcase_08 AC 11 ms
4,380 KB
testcase_09 AC 16 ms
4,380 KB
権限があれば一括ダウンロードができます

ソースコード

diff # 

#include <stdio.h>
#include <stdint.h>
#include <math.h>

#pragma region type
/* signed integer */
typedef   int8_t      i8;
typedef   int16_t     i16;
typedef   int32_t     i32;
typedef   int64_t     i64;
typedef __int128_t    i128;
/* unsigned integer */
typedef   uint8_t     u8;
typedef   uint16_t    u16;
typedef   uint32_t    u32;
typedef   uint64_t    u64;
typedef __uint128_t   u128;
/* floating point number */
typedef   float       f32;
typedef   double      f64;
typedef   long double f80;
#pragma endregion type


#pragma region m64
typedef uint64_t m64;
m64 _one_m64(u64 mod) { return (u64)-1ull % mod + 1; }
m64 _r2_m64(u64 mod) { return (u128)(i128)-1 % mod + 1; }
m64 _inv_m64(u64 mod) {
  m64 inv = mod;
  for (int i = 0; i < 5; i++) inv *= 2 - inv * mod;
  return inv;
}
m64 _reduce_m64(u128 a, m64 inv, u64 mod) {
  u64 y = (u64)(a >> 64) - (u64)(((u128)((u64)a * inv) * mod) >> 64);
  return (i64)y < 0 ? y + mod : y;
}
m64 to_m64(u64 a, m64 r2, m64 inv, u64 mod) { return _reduce_m64((u128)a * r2, inv, mod); }
m64 add_m64(m64 A, m64 B, u64 mod) {
  return A + B >= mod ? A + B - mod:  A + B;
}
m64 sub_m64(m64 A, m64 B, u64 mod) {
  return A >= B ? A - B : mod + A - B;
}
m64 mul_m64(m64 A, m64 B, m64 inv, u64 mod) { return _reduce_m64((u128)A * B, inv, mod); }
#pragma endregion m64

int jacobi(i64 a, u64 n)
{
  u64 t;
  int j = 1;
  while (a)
  {
    if (a < 0)
    {
      a = -a;
      if ((n & 3) == 3)
        j = -j;
    }
    int ba = __builtin_ctzll(a);
    a >>= ba;
    if (((n & 7) == 3 || (n & 7) == 5) && (ba & 1))
      j = -j;
    if ((a & n & 3) == 3)
      j = -j;
    t = a;
    a = n;
    n = t;
    a %= n;
    if (a > n / 2)
      a -= n;
  }
  return n == 1 ? j : 0;
}

int is_prime(const u64 n)
{
  if (n <= 1)
    return 0;
  if (n <= 3)
    return 1;
  if (!(n & 1))
    return 0;

  const m64 one = -1ULL % n + 1;
  const m64 r2 = (u128)(i128)-1 % n + 1;
  const m64 m = _inv_m64(n);

  {
    uint64_t d = (n - 1) << __builtin_clzll(n - 1);
    m64 t = one << 1;
    if (t >= n)
      t -= n;
    for (d <<= 1; d; d <<= 1)
    {
      t = mul_m64(t, t, m, n);
      if (d >> 63)
      {
        t <<= 1;
        if (t >= n)
          t -= n;
      }
    }
    if (t != one)
    {
      u64 x = (n - 1) & -(n - 1);
      m64 mone = n - one;
      for (x >>= 1; t != mone; x >>= 1)
      {
        if (x == 0)
          return 0;
        t = mul_m64(t, t, m, n);
      }
    }
  }

  {
    i64 D = 5;
    int i;
    for (i = 0; jacobi(D, n) != -1 && i < 64; i++)
    {
      if (i == 32)
      {
        u32 k = round(sqrtl(n));
        if (k * k == n)
          return 0;
      }
      if (i & 1)
        D -= 2;
      else
        D += 2;
      D = -D;
    }

    m64 Q = to_m64(D < 0 ? (1 - D) / 4 % n : n - (D - 1) / 4 % n, r2, m, n);

    m64 u, v, Qn;
    u64 k = (n + 1) << __builtin_clzll(n + 1);
    u = one;
    v = one;
    Qn = Q;
    D %= (i64)n;
    D = to_m64(D < 0 ? n + D : D, r2, m, n);
    for (k <<= 1; k; k <<= 1)
    {
      u = mul_m64(u, v, m, n);
      v = sub_m64(mul_m64(v, v, m, n), add_m64(Qn, Qn, n), n);
      Qn = mul_m64(Qn, Qn, m, n);
      if (k >> 63)
      {
        u64 uu = add_m64(u, v, n);
        if (uu & 1)
          uu += n;
        uu >>= 1;
        v = add_m64(mul_m64(D, u, m, n), v, n);
        if (v & 1)
          v += n;
        v >>= 1;
        u = uu;
        Qn = mul_m64(Qn, Q, m, n);
      }
    }

    if (u == 0 || v == 0)
      return 1;
    u64 x = (n + 1) & ~n;
    for (x >>= 1; x; x >>= 1)
    {
      u = mul_m64(u, v, m, n);
      v = sub_m64(mul_m64(v, v, m, n), add_m64(Qn, Qn, n), n);
      if (v == 0)
        return 1;
      Qn = mul_m64(Qn, Qn, m, n);
    }
  }

  return 0;
}

int main()
{
  int T; scanf("%d", &T);
  while (T--) {
    unsigned long long int x;
    scanf("%llu", &x);
    printf("%llu %d\n", x, is_prime(x));
  }
  return 0;
}
0