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

Nunciature

Why not use constexpr? Take a look here : cpptruths.blogspot.no/2011/07/… – Beta

So, you're being trolled by a C++ guru and come to SO with the ... homework... :) – Epiphysis

@Beta I dunno, but not all compilers that claim to support C++11 have constexpr. (I'm looking at you VS2012...) – Howling

@olevegard, heh, I just forgot about its existance. Could you post it as answer? – Nunciature

@sehe, believe you or not, I did wrote this code myself and It was my own idea. Am I trolled by me? Maybe. But I wouldn't say he(me) is a guru – Nunciature

You could cut your depth by half fairly easily if you special case 2 and only check for odd factors in the recursion. – Inconsecutive

@Casey: yes it's good idea, but unfortunately it will be not ehough. – Nunciature

@Nunciature Lol. That counts :) – Epiphysis

@Nunciature Out of curiosity; why can't you change the compilation string? Personal challenge? – Beta

@olevegard, I want to prepare some code that could be used in online judges. – Nunciature

Couldn't you include a look-up table? Probably not for the whole thing, but for numbers up to some N. (E.g. as partial specializations of some template for prime numbers) – Declivity

@DyP, I'm not sure how can I use that to make it help. I already can calculate everything approx. up to 10^6, but It will extremely increase source code to print all primes up to 10^6, – Nunciature

One obvious way to reduce the template depth is to test for multiple of 2 separately, and then only test divisibility by odd numbers. You can even go further by eliminating e.g. multiples of 3 and 5, although then the updating gets more complex. – Coverdale

Too bad, using the multiples of 2 technique only helps until 10**9+6 : coliru.stacked-crooked.com/…. – Erivan

take a look at this question: https://mcmap.net/q/740634/-what-39-s-the-most-efficient-tail-recursive-prime-verification-function-known/1000282 – Dealings

Note that constexpr is not a guarantee of compile-time evaluation. It may or may not be evaluated at compile time. Unless it is used in a place where the value is required at compile-time (i.e. template argument). – Decry

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.

Inconsecutive

"Apologies for the horrible formating..." don't worry, is thousands of times better than the equivalent template-meta-programming format :) – Ottillia

I like that. I use it checking even numbers too and it looks a bit more cleaner. Thanks. – Nunciature

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

Mummer

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

Vaunting

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

Beta

Hmm, I've noticed its depth is limited to (even limit is lower): ideone.com/VGJxnr – Nunciature

Read the article, you can set the limit with -fconstexpr-depth=9300 If you set the depth manually, constexpr can "go deeper" than template meta programming. At least according to the article, haven't tried myself. – Beta

I can't change compilation string. If I could I would change it for my own code. – Nunciature

coliru.stacked-crooked.com It worked using this compiler for me : coliru.stacked-crooked.com/… – Beta

It works because 7 is quite low number. It doesn't compiled for 1000000007 for instance – Nunciature

@RiaD, it works, just increase the -fconstexpr-depth=930000 – Beta

@olevegard: "I can't change compilation string. If I could I would change it for my own code." - "I can't increase the depth limit." – Prolongation

@Prolongation Yea, I kinda forgot about that. But it is worth mentioning to others who might see this in the future. – Beta

One thing: you can lower the number of calls to functions if you put a special case for 2 in is_prime_func and have your recursive use c+2 instead of c+1. – Erivan

You can't add compilation flags? Why not? – Kramer

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.

Holmen

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

Voletta

Well, it basically the same as in the OP and doesn't solve desired problem ideone.com/u45RNt – Nunciature

Enums are always evaluated in compile time so as IsPrime<17>::Value. Are you looking for something different ? – Voletta

I need it to work for big numbers without hitting limit of instantiation depth. – Nunciature