Prime Number Utilities
(weilycoder/number-theory/prime.hpp)

• View this file on GitHub
• Last update: 2025-11-07 09:14:02+08:00
• Include: #include "weilycoder/number-theory/prime.hpp"

Depends on

• Modular Arithmetic Utilities (weilycoder/number-theory/mod_utility.hpp)

Required by

• Functions for factorizing numbers using Pollard's Rho algorithm (weilycoder/number-theory/factorize.hpp)
• Functions to find primitive roots modulo a prime (weilycoder/number-theory/primitive_root.hpp)
• Number Theoretic Transform (NTT) (weilycoder/poly/ntt.hpp)
• weilycoder/poly/ntt_convolve.hpp

Verified with

• test/convolution_mod.test.cpp
• test/factorize.test.cpp
• test/primality_test.test.cpp
• test/primitive_root.test.cpp

Code

#ifndef WEILYCODER_PRIME_HPP
#define WEILYCODER_PRIME_HPP

/**
 * @file prime.hpp
 * @brief Prime Number Utilities
 */

#include "mod_utility.hpp"
#include <cstdint>
#include <type_traits>

namespace weilycoder {
/**
 * @brief Miller-Rabin primality test for a given base.
 * @tparam bit32 If true, won't use 128-bit arithmetic. You should ensure that
 *         all inputs are small enough to avoid overflow (i.e. bit-32).
 * @tparam base The base to test.
 * @param n The number to test for primality.
 * @param d An odd component of n-1 (n-1 = d * 2^s).
 * @param s The exponent of 2 in the factorization of n-1.
 * @return true if n is probably prime for the given base, false if composite.
 */
template <bool bit32, uint64_t base>
constexpr bool miller_rabin_test(uint64_t n, uint64_t d, uint32_t s) {
  uint64_t x = mod_pow<bit32>(base, d, n);
  if (x == 0 || x == 1 || x == n - 1)
    return true;
  for (uint32_t r = 1; r < s; ++r) {
    x = mod_mul<bit32>(x, x, n);
    if (x == n - 1)
      return true;
  }
  return false;
}

/**
 * @brief Variadic template to test multiple bases in Miller-Rabin test.
 * @tparam bit32 If true, won't use 128-bit arithmetic. You should ensure that
 *         all inputs are small enough to avoid overflow (i.e. bit-32).
 * @tparam base The first base to test.
 * @tparam Rest The remaining bases to test.
 * @param n The number to test for primality.
 * @param d An odd component of n-1 (n-1 = d * 2^s).
 * @param s The exponent of 2 in the factorization of n-1.
 * @return true if n is probably prime for all given bases, false if composite.
 */
template <bool bit32, uint64_t base, uint64_t... Rest>
constexpr std::enable_if_t<(sizeof...(Rest) != 0), bool>
miller_rabin_test(uint64_t n, uint64_t d, uint32_t s) {
  return miller_rabin_test<bit32, base>(n, d, s) &&
         miller_rabin_test<bit32, Rest...>(n, d, s);
}

/**
 * @brief Miller-Rabin primality test using multiple bases.
 * @tparam bit32 If true, won't use 128-bit arithmetic. You should ensure that
 *         all inputs are small enough to avoid overflow (i.e. bit-32).
 * @tparam bases The bases to test.
 * @param n The number to test for primality.
 * @return true if n is probably prime, false if composite.
 */
template <bool bit32, uint64_t... bases> constexpr bool miller_rabin(uint64_t n) {
  if (n < 2)
    return false;
  if (n == 2 || n == 3)
    return true;
  if (n % 2 == 0)
    return false;

  uint64_t d = n - 1, s = 0;
  for (; d % 2 == 0; d /= 2)
    ++s;

  return miller_rabin_test<bit32, bases...>(n, d, s);
}

/**
 * @brief Miller-Rabin primality test optimized for 64-bit integers.
 *        Uses a fixed set of bases that guarantee correctness
 *        for 64-bit integers.
 * @param n The number to test for primality.
 * @return true if n is prime, false if not prime.
 */
constexpr bool miller_rabin64(uint64_t n) {
  return miller_rabin<false, 2, 325, 9375, 28178, 450775, 9780504, 1795265022>(n);
}

/**
 * @brief Miller-Rabin primality test optimized for 32-bit integers.
 *        Uses a fixed set of bases that guarantee correctness
 *        for 32-bit integers.
 * @param n The number to test for primality.
 * @return true if n is prime, false if not prime.
 */
constexpr bool miller_rabin32(uint32_t n) { return miller_rabin<true, 2, 7, 61>(n); }

constexpr bool is_prime(uint64_t n) {
  if (n <= UINT32_MAX)
    return miller_rabin32(static_cast<uint32_t>(n));
  return miller_rabin64(n);
}

/**
 * @brief Compile-time primality test for a given integer.
 * @tparam x The integer to test for primality.
 * @return true if x is prime, false otherwise.
 */
template <uint64_t x> constexpr bool is_prime() { return is_prime(x); }
} // namespace weilycoder

#endif

#line 1 "weilycoder/number-theory/prime.hpp"



/**
 * @file prime.hpp
 * @brief Prime Number Utilities
 */

#line 1 "weilycoder/number-theory/mod_utility.hpp"



/**
 * @file mod_utility.hpp
 * @brief Modular Arithmetic Utilities
 */

#include <cstdint>

namespace weilycoder {
using u128 = unsigned __int128;

template <bool bit32 = false>
constexpr uint64_t mod_add(uint64_t a, uint64_t b, uint64_t modulus) {
  if constexpr (bit32) {
    uint64_t res = a + b;
    if (res >= modulus)
      res -= modulus;
    return res;
  } else {
    u128 res = static_cast<u128>(a) + b;
    if (res >= modulus)
      res -= modulus;
    return res;
  }
}

template <bool bit32 = false>
constexpr uint64_t mod_sub(uint64_t a, uint64_t b, uint64_t modulus) {
  if constexpr (bit32) {
    uint64_t res = (a >= b) ? (a - b) : (modulus + a - b);
    return res;
  } else {
    u128 res = (a >= b) ? (a - b) : (static_cast<u128>(modulus) + a - b);
    return res;
  }
}

/**
 * @brief Perform modular multiplication for 64-bit integers.
 * @tparam bit32 If true, won't use 128-bit arithmetic. You should ensure that
 *         all inputs are small enough to avoid overflow (i.e. bit-32).
 * @param a The first multiplicand.
 * @param b The second multiplicand.
 * @param modulus The modulus.
 * @return (a * b) % modulus
 */
template <bool bit32 = false>
constexpr uint64_t mod_mul(uint64_t a, uint64_t b, uint64_t modulus) {
  if constexpr (bit32)
    return a * b % modulus;
  else
    return static_cast<u128>(a) * b % modulus;
}

/**
 * @brief Perform modular multiplication for 64-bit integers with a compile-time
 *        modulus.
 * @tparam Modulus The modulus.
 * @tparam bit32 If true, won't use 128-bit arithmetic. You should ensure that
 *         all inputs are small enough to avoid overflow (i.e. bit-32).
 * @param a The first multiplicand.
 * @param b The second multiplicand.
 * @return (a * b) % Modulus
 */
template <uint64_t Modulus, bool bit32 = false>
constexpr uint64_t mod_mul(uint64_t a, uint64_t b) {
  if constexpr (bit32)
    return a * b % Modulus;
  else
    return static_cast<u128>(a) * b % Modulus;
}

/**
 * @brief Perform modular exponentiation for 64-bit integers.
 * @tparam bit32 If true, won't use 128-bit arithmetic. You should ensure that
 *         all inputs are small enough to avoid overflow (i.e. bit-32).
 * @param base The base number.
 * @param exponent The exponent.
 * @param modulus The modulus.
 * @return (base^exponent) % modulus
 */
template <bool bit32 = false>
constexpr uint64_t mod_pow(uint64_t base, uint64_t exponent, uint64_t modulus) {
  uint64_t result = 1 % modulus;
  base %= modulus;
  while (exponent > 0) {
    if (exponent & 1)
      result = mod_mul<bit32>(result, base, modulus);
    base = mod_mul<bit32>(base, base, modulus);
    exponent >>= 1;
  }
  return result;
}
} // namespace weilycoder


#line 11 "weilycoder/number-theory/prime.hpp"
#include <type_traits>

namespace weilycoder {
/**
 * @brief Miller-Rabin primality test for a given base.
 * @tparam bit32 If true, won't use 128-bit arithmetic. You should ensure that
 *         all inputs are small enough to avoid overflow (i.e. bit-32).
 * @tparam base The base to test.
 * @param n The number to test for primality.
 * @param d An odd component of n-1 (n-1 = d * 2^s).
 * @param s The exponent of 2 in the factorization of n-1.
 * @return true if n is probably prime for the given base, false if composite.
 */
template <bool bit32, uint64_t base>
constexpr bool miller_rabin_test(uint64_t n, uint64_t d, uint32_t s) {
  uint64_t x = mod_pow<bit32>(base, d, n);
  if (x == 0 || x == 1 || x == n - 1)
    return true;
  for (uint32_t r = 1; r < s; ++r) {
    x = mod_mul<bit32>(x, x, n);
    if (x == n - 1)
      return true;
  }
  return false;
}

/**
 * @brief Variadic template to test multiple bases in Miller-Rabin test.
 * @tparam bit32 If true, won't use 128-bit arithmetic. You should ensure that
 *         all inputs are small enough to avoid overflow (i.e. bit-32).
 * @tparam base The first base to test.
 * @tparam Rest The remaining bases to test.
 * @param n The number to test for primality.
 * @param d An odd component of n-1 (n-1 = d * 2^s).
 * @param s The exponent of 2 in the factorization of n-1.
 * @return true if n is probably prime for all given bases, false if composite.
 */
template <bool bit32, uint64_t base, uint64_t... Rest>
constexpr std::enable_if_t<(sizeof...(Rest) != 0), bool>
miller_rabin_test(uint64_t n, uint64_t d, uint32_t s) {
  return miller_rabin_test<bit32, base>(n, d, s) &&
         miller_rabin_test<bit32, Rest...>(n, d, s);
}

/**
 * @brief Miller-Rabin primality test using multiple bases.
 * @tparam bit32 If true, won't use 128-bit arithmetic. You should ensure that
 *         all inputs are small enough to avoid overflow (i.e. bit-32).
 * @tparam bases The bases to test.
 * @param n The number to test for primality.
 * @return true if n is probably prime, false if composite.
 */
template <bool bit32, uint64_t... bases> constexpr bool miller_rabin(uint64_t n) {
  if (n < 2)
    return false;
  if (n == 2 || n == 3)
    return true;
  if (n % 2 == 0)
    return false;

  uint64_t d = n - 1, s = 0;
  for (; d % 2 == 0; d /= 2)
    ++s;

  return miller_rabin_test<bit32, bases...>(n, d, s);
}

/**
 * @brief Miller-Rabin primality test optimized for 64-bit integers.
 *        Uses a fixed set of bases that guarantee correctness
 *        for 64-bit integers.
 * @param n The number to test for primality.
 * @return true if n is prime, false if not prime.
 */
constexpr bool miller_rabin64(uint64_t n) {
  return miller_rabin<false, 2, 325, 9375, 28178, 450775, 9780504, 1795265022>(n);
}

/**
 * @brief Miller-Rabin primality test optimized for 32-bit integers.
 *        Uses a fixed set of bases that guarantee correctness
 *        for 32-bit integers.
 * @param n The number to test for primality.
 * @return true if n is prime, false if not prime.
 */
constexpr bool miller_rabin32(uint32_t n) { return miller_rabin<true, 2, 7, 61>(n); }

constexpr bool is_prime(uint64_t n) {
  if (n <= UINT32_MAX)
    return miller_rabin32(static_cast<uint32_t>(n));
  return miller_rabin64(n);
}

/**
 * @brief Compile-time primality test for a given integer.
 * @tparam x The integer to test for primality.
 * @return true if x is prime, false otherwise.
 */
template <uint64_t x> constexpr bool is_prime() { return is_prime(x); }
} // namespace weilycoder

Back to top page