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 { };

Answer 4

You can take a look at constexpr. It has a lot friendlier syntax than template meta programming ( at least if you're unfamiliar with templates like me. ) You can't use if's, or any loops. But with recursion and the tenary opeartor you can do pretty much everything you can with template meta programming, and it usually runs faster too.

http://cpptruths.blogspot.no/2011/07/want-speed-use-constexpr-meta.html

Here's a working example using an online compiler : http://coliru.stacked-crooked.com/view?id=6bc10e71b8606dd2980c0c5dd982a3c0-6fbdb8a7476ab90c2bd2503cd4005881

Since it is executed at compile time, you can do a static assert to test it.

static_assert(is_prime_func(x), "...");

The assert will fail if x is not a prime, meaning compilation will fail. If x is a prime, the compilation will succeed but no output will be generated.

If you want to check really large numbers you can increase the constexpr depth

-fconstexpr-depth=930000

I haven't tested how large numbers it supports, but I assume it varies from compiler to compiler.

If you want to test it out for yourself :

#include <cstdio>

constexpr bool is_prime_recursive(size_t number, size_t c)
{
  return (c*c > number) ? true : 
           (number % c == 0) ? false : 
              is_prime_recursive(number, c+1);
}

constexpr bool is_prime_func(size_t number)
{
  return (number <= 1) ? false : is_prime_recursive(number, 2);
}

int main(void)
{
   static_assert(is_prime_func(7), "...");  // Computed at compile-time
}

Compiling

g++ -std=c++11 -O2 -Wall -pedantic -pthread main.cpp -std=c++11 -fconstexpr-depth=9300 && ./a.out

Answer 5

From a language point of view, the solution is to increase the depth limit. The program is correct except it requires "too many" iterations. But you have stated you don't want to increase it. (It appears that the required template depth is (sqrt(N) + C) where C is a small constant, so for a max depth of 900 on your system your current implementation will work up to 810000.)

I can think of two strategies to increase the upper range limit:

1. Improve the algorithm. If you only check odd factors, you cut the number of iterations in half. The upper limit goes up by a factor of four. This still isn't anywhere near a billion, but you can certainly get there by more closely approximating the ideal sieve.

2. Use a typedef declaration to pre-evaluate part of the metafunction, and rely on your compiler's memoization policy to stop that part from being reevaluated at full depth.

   This strategy won't work for metafunctions that rely heavily on results from previous iterations, but in your case you can check the last 900 factors, and then a check for the last 1800 factors will automatically use a cached copy of the result from the last 900. This isn't specified by the C++ standard, and strictly isn't portable, but on the other hand neither is anything regarding these recursion limits.

Answer 6

C++ without constexpr, IsPrime::Value give out the compile-time result. The trick is to iteratively try to divide by i=3,5,7,... until i*i>n

template <int n, int i, int b> struct IsPrimeIter;

template <int n, int i>
struct IsPrimeIter<n, i, 0> {
    enum _ { Value = 0 };
};

template <int n, int i>
struct IsPrimeIter<n, i, 1> {
    enum _ { Value = 1 };
};

template <int n, int i>
struct IsPrimeIter<n, i, 2> {
    enum _ { Value = IsPrimeIter<n, i+2, 
                             (i * i > n) ? 1 : 
                             n % i == 0 ? 0 : 2>::Value };
};

template <int n>
struct IsPrime {
    enum _ { Value = n <= 1 ? false:
                     (n == 2 || n == 3) ? true:
                     (n % 2 == 0) ? false :
                     IsPrimeIter<n, 3, 2>::Value };
};