Miller-Rabin 素数判定
(number-theory/miller-rabin.hpp)

• View this file on GitHub
• Last update: 2026-07-25 02:01:37+09:00
• Include: #include "number-theory/miller-rabin.hpp"

64 bit 整数の素数判定を行う．

• MillerRabin::is_prime(N)：$N$ が素数であるか判定する．$N\leq 10^{18}$ とし，$N<2$ のとき false を返す．

64 bit 整数に対して決定的となる基底を用いる．時間計算量は $O(\log N)$ 回の剰余乗算．

Required by

• Pollard's rho algorithm (number-theory/pollard-rho.hpp)
• 原始根 (number-theory/primitive-root.hpp)

Verified with

• verify/number-theory/LC_factorize.test.cpp
• verify/number-theory/LC_primality_test.test.cpp
• verify/number-theory/LC_primitive_root.test.cpp
• verify/number-theory/UNIT_pollard_rho_divisors.test.cpp

Code

#pragma once

namespace MillerRabin {
using u64 = uint64_t;
using u128 = __uint128_t;

namespace internal {
u64 multiply_mod(u64 a, u64 b, u64 mod) { return u128(a) * b % mod; }

u64 power_mod(u64 a, u64 n, u64 mod) {
  u64 ret = 1;
  while (n) {
    if (n & 1) ret = multiply_mod(ret, a, mod);
    a = multiply_mod(a, a, mod);
    n >>= 1;
  }
  return ret;
}
};  // namespace internal

bool is_prime(long long n) {
  if (n < 2) return false;
  u64 x = n;
  for (u64 p : {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37}) {
    if (x % p == 0) return x == p;
  }

  int s = __builtin_ctzll(x - 1);
  u64 d = (x - 1) >> s;
  for (u64 a : {2, 325, 9375, 28178, 450775, 9780504, 1795265022}) {
    if (a % x == 0) continue;
    u64 y = internal::power_mod(a % x, d, x);
    if (y == 1 || y == x - 1) continue;
    bool composite = true;
    for (int r = 1; r < s; r++) {
      y = internal::multiply_mod(y, y, x);
      if (y == x - 1) {
        composite = false;
        break;
      }
    }
    if (composite) return false;
  }
  return true;
}
};  // namespace MillerRabin

/**
 * @brief Miller-Rabin 素数判定
 * @docs docs/number-theory/miller-rabin.md
 */

#line 2 "number-theory/miller-rabin.hpp"

namespace MillerRabin {
using u64 = uint64_t;
using u128 = __uint128_t;

namespace internal {
u64 multiply_mod(u64 a, u64 b, u64 mod) { return u128(a) * b % mod; }

u64 power_mod(u64 a, u64 n, u64 mod) {
  u64 ret = 1;
  while (n) {
    if (n & 1) ret = multiply_mod(ret, a, mod);
    a = multiply_mod(a, a, mod);
    n >>= 1;
  }
  return ret;
}
};  // namespace internal

bool is_prime(long long n) {
  if (n < 2) return false;
  u64 x = n;
  for (u64 p : {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37}) {
    if (x % p == 0) return x == p;
  }

  int s = __builtin_ctzll(x - 1);
  u64 d = (x - 1) >> s;
  for (u64 a : {2, 325, 9375, 28178, 450775, 9780504, 1795265022}) {
    if (a % x == 0) continue;
    u64 y = internal::power_mod(a % x, d, x);
    if (y == 1 || y == x - 1) continue;
    bool composite = true;
    for (int r = 1; r < s; r++) {
      y = internal::multiply_mod(y, y, x);
      if (y == x - 1) {
        composite = false;
        break;
      }
    }
    if (composite) return false;
  }
  return true;
}
};  // namespace MillerRabin

/**
 * @brief Miller-Rabin 素数判定
 * @docs docs/number-theory/miller-rabin.md
 */

Back to top page