Compile time prime checking

Question

I need to check is some integer prime in compile time (to put the boolean value as template argument).

I've write code that do it well:

#include <type_traits>
namespace impl {
    template <int n, long long i>
    struct PrimeChecker {
        typedef typename std::conditional<
                    (i * i > n),
                    std::true_type,
                    typename std::conditional<
                        n % i == 0,
                        std::false_type,
                        typename PrimeChecker<n, (i * i > n ) ? -1 : i + 1>::type
                    >::type
                >::type type;
    };
    template <int n>
    struct PrimeChecker<n, -1> {
        typedef void type;
    };
} // namespace impl
template<int n>
struct IsPrime {
    typedef typename impl::PrimeChecker<n, 2>::type type;
};

template<>
struct IsPrime<1> : public std::false_type {
};

It works for numbers to ~1000000 and fails with error for 10⁹

prog.cpp:15:23: error: template instantiation depth exceeds maximum of 900 (use -ftemplate-depth= to increase the maximum) instantiating ‘struct impl::PrimeChecker<1000000000, 901ll>’
               >::type type;
                       ^
prog.cpp:15:23:   recursively required from ‘struct impl::PrimeChecker<1000000000, 3ll>’
prog.cpp:15:23:   required from ‘struct impl::PrimeChecker<1000000000, 2ll>’
prog.cpp:24:54:   required from ‘struct IsPrime<1000000000>’
prog.cpp:32:41:   required from here

I can't increase the depth limit. Is it somehow possible to decrease depth I use?

Thing I want to achive: I need to check is constant prime in compile time without changing compilation string with template depth limit 900 and constexpr depth limit 512. (default for my g++). It should work for all positive int32's or at least for numbers up to 10⁹+9

Answer 1

You can change the space requirement from linear to logarithmic by splitting the range by halves using a divide-and-conquer algorithm. This method uses divide-and-conquer, and only tests odd factors (Live at Coliru):

namespace detail {
using std::size_t;

constexpr size_t mid(size_t low, size_t high) {
  return (low + high) / 2;
}

// precondition: low*low <= n, high*high > n.
constexpr size_t ceilsqrt(size_t n, size_t low, size_t high) {
  return low + 1 >= high
    ? high
    : (mid(low, high) * mid(low, high) == n)
      ? mid(low, high)
      : (mid(low, high) * mid(low, high) <  n)
        ? ceilsqrt(n, mid(low, high), high)
        : ceilsqrt(n, low, mid(low, high));
}

// returns ceiling(sqrt(n))
constexpr size_t ceilsqrt(size_t n) {
  return n < 3
    ? n
    : ceilsqrt(n, 1, size_t(1) << (std::numeric_limits<size_t>::digits / 2));
}


// returns true if n is divisible by an odd integer in
// [2 * low + 1, 2 * high + 1).
constexpr bool find_factor(size_t n, size_t low, size_t high)
{
  return low + 1 >= high
    ? (n % (2 * low + 1)) == 0
    : (find_factor(n, low, mid(low, high))
       || find_factor(n, mid(low, high), high));
}
}

constexpr bool is_prime(std::size_t n)
{
  using detail::find_factor;
  using detail::ceilsqrt;

  return n > 1
    && (n == 2
        || (n % 2 == 1
            && (n == 3
                || !find_factor(n, 1, (ceilsqrt(n) + 1) / 2))));
}

EDIT: Use compile-time sqrt to bound search space to ceiling(sqrt(n)), instead of n / 2. Now can compute is_prime(100000007) as desired (and is_prime(1000000000039ULL) for that matter) on Coliru without exploding.

Apologies for the horrible formatting, I still haven't found a comfortable style for C++11's tortured constexpr sub-language.

EDIT: Cleanup code: replace macro with another function, move implementation detail into a detail namespace, steal indentation style from Pablo's answer.

Answer 2

Here's my attempt at this. Using constexpr and the deterministic variant of the Miller-Rabin primality test for numbers up to 4,759,123,141 (which should cover all uint32s, but you can easily change the primer-checker set to cover a larger range)

#include <cstdint>

constexpr uint64_t ct_mod_sqr(uint64_t a, uint64_t m)
{
  return (a * a) % m;
}

constexpr uint64_t ct_mod_pow(uint64_t a, uint64_t n, uint64_t m)
{
  return (n == 0)
    ? 1
    : (ct_mod_sqr(ct_mod_pow(a, n/2, m), m) * ((n & 1) ? a : 1)) % m;
}

constexpr bool pass_prime_check_impl(uint64_t x, uint32_t n1, uint32_t s1)
{
  return (s1 == 0)
    ? false
    : (x == 1)
      ? false
      : (x == n1)
        ? true
        : pass_prime_check_impl(ct_mod_sqr(x, n1 + 1), n1, s1 - 1)
    ;
}

constexpr bool pass_prime_check_impl(uint32_t a, uint32_t n1, uint32_t s1, uint32_t d, uint64_t x)
{
  return (x == 1) || (x == n1)
    ? true
    : pass_prime_check_impl(ct_mod_sqr(x, n1 + 1), n1, s1)
    ;
}

constexpr bool pass_prime_check_impl(uint32_t a, uint32_t n1, uint32_t s1, uint32_t d)
{
  return pass_prime_check_impl(a, n1, s1, d,
                               ct_mod_pow(a, d, n1 + 1));
}

constexpr bool pass_prime_check_impl(uint32_t n, uint32_t a)
{
  return pass_prime_check_impl(a, n - 1,
                               __builtin_ctz(n - 1) - 1,
                               (n - 1) >> __builtin_ctz(n - 1));
}

constexpr bool pass_prime_check(uint32_t n, uint32_t p)
{
  return (n == p)
    ? true
    : pass_prime_check_impl(n, p);
}

constexpr bool is_prime(uint32_t n)
{
  return (n == 2)
    ? true
    : (n % 2 == 0)
      ? false
      : (pass_prime_check(n, 2) &&
         pass_prime_check(n, 7) &&
         pass_prime_check(n, 61))
    ;
}

int main()
{
  static_assert(is_prime(100000007),  "100000007 is a prime!");
  static_assert(is_prime(1000000007), "1000000007 is a prime!");
  static_assert(is_prime(1000000009), "1000000009 is a prime!");
  static_assert(!is_prime(1000000011), "1000000011 is not a prime!");
  return 0;
}

Answer 3

constexpr are probably easier to deal with, but there's no real problem doing this with pure template instantiation.

UPDATE: Fixed Newton-Raphson integer square root

This code is suboptimal -- dropping all the test divisions by even numbers (and maybe even by multiples of three) would obviously speed up compile times -- but it works and even with a prime about 10¹⁰ gcc uses less than 1GB of RAM.

#include <type_traits>

template<typename a, typename b> struct both
  : std::integral_constant<bool, a::value && b::value> { };

template<long long low, long long spread, long long n>
struct HasNoFactor
  : both<typename HasNoFactor<low, spread/2, n>::type,
         typename HasNoFactor<low+spread/2, (spread + 1)/2, n>::type> { };

template<long long low, long long n>
struct HasNoFactor<low, 0, n> : std::true_type { };

template<long long low, long long n>
struct HasNoFactor<low, 1, n>
  : std::integral_constant<bool, n % low != 0> { };

// Newton-Raphson computation of floor(sqrt(n))

template<bool done, long long n, long long g>
struct ISqrtStep;

template<long long n, long long g = n, long long h = (n + 1) / 2, bool done = (g <= h)>
struct ISqrt;

template<long long n, long long g, long long h>
struct ISqrt<n, g, h, true> : std::integral_constant<long long, g> { };

template<long long n, long long g, long long h>
struct ISqrt<n, g, h, false> : ISqrt<n, h, (h + n / h) / 2> { };

template<long long n>
struct IsPrime : HasNoFactor<2, ISqrt<n>::value - 1, n> { };

template<> struct IsPrime<0> : std::false_type { };
template<> struct IsPrime<1> : std::false_type { };