constexpr is_prime using Miller test

is_prime.cpp

#include <cstdint>

using u64 = std::uint64_t;

// a + b
constexpr u64 mod_add(u64 a, u64 b, u64 m) {
  if (~a < b) {
    return (a + b) % m + ~m + 1;
  } else {
    return a + b;
  }
}
// a * 2^i
constexpr u64 mod_power(u64 a, int i, u64 m) {
  if (i < 1) {
    return a;
  } else {
    return mod_power(mod_add(a, a, m), i - 1, m);
  }
}
constexpr u64 hiword(u64 x) {
  return x >> 32;
}
constexpr u64 loword(u64 x) {
  return x & hiword(~0);
}
// a * b
constexpr void mod_mul(u64& a, u64 b, u64 m) {
  u64 x = 1;
  if (m < (x << 32)) {
    a = (a % m) * (b % m);
  } else {
    x = 0;
    x = mod_add(x, hiword(a) * hiword(b), m);
    x = mod_power(x, 32, m);
    x = mod_add(x, loword(a) * hiword(b), m);
    x = mod_add(x, hiword(a) * loword(b), m);
    x = mod_power(x, 32, m);
    a = mod_add(x, loword(a) * loword(b), m);
  }
}
constexpr bool even(u64 n) {
  return n % 2 == 0;
}
constexpr bool miller_test(u64 a, u64 s, u64 d, u64 n) {
  u64 x = 1;
  while (d != 0) {
    if (even(d)) {
      mod_mul(a, a, n);
      d /= 2;
    } else {
      mod_mul(x, a, n);
      d -= 1;
    }
  }
  // x = a^d mod n
  if (x % n == 1) {
    return true;
  }
  for (u64 r = 0; r < s; ++r) {
    // x = a^(d * 2^r) mod n
    if (x % n == n - 1) {
      return true;
    }
    mod_mul(x, x, n);
  }
  return false;
}
constexpr bool is_prime(u64 n) {
  if (n < 2 || even(n)) {
    return n == 2;
  }
  u64 s = 0, d = n - 1;
  while (even(d)) {
    d /= 2;
    ++s;
  }
  // n - 1 = d * 2^s
  u64 x[] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37};
  for (u64 a : x) {
    if (a < n && !miller_test(a, s, d, n)) {
      return false;
    }
  }
  return true;
}