Compute nth prime at compile time [closed]

Question

This question is unlikely to help any future visitors; it is only relevant to a small geographic area, a specific moment in time, or an extraordinarily narrow situation that is not generally applicable to the worldwide audience of the internet. For help making this question more broadly applicable, visit the help center.

Closed 7 years ago.

The C++11 features, with constexpr and template argument packs, should in my opinion be strong enough to perform some rather complex computations. One possible example for which I have a practical application is the computation of the nth prime at compile time.

I'm asking for ways to implement this computation. If more than one solution are proposed, it might be interesting to compare them.

To give you an idea of my performance expectations: I'd hope for some code which can find the 512th prime (which is 3671) in less than one second compile time on reasonable desktop hardware.

Answer 1

I implemented the simplest possible way, not using templates at all and it works:

constexpr bool isPrimeLoop(int i, int k) {
    return (k*k > i)?true:(i%k == 0)?false:isPrimeLoop(i, k + 1);
}

constexpr bool isPrime(int i) {
    return isPrimeLoop(i, 2);
}

constexpr int nextPrime(int k) {
    return isPrime(k)?k:nextPrime(k + 1);
}

constexpr int getPrimeLoop(int i, int k) {
// i - nr of primes to advance
// k - some starting prime
    return (i == 0)?k:getPrimeLoop(i - 1, nextPrime(k + 1));
}

constexpr int getPrime(int i) {
    return getPrimeLoop(i, 2);
}

static_assert(getPrime(511) == 3671, "computed incorrectly");

It needs increased constexpr-depth a bit, but it fits in time easily:

$ time g++ -c -std=c++11 vec.cpp -fconstexpr-depth=600

real    0m0.093s
user    0m0.080s
sys 0m0.008s

The following trick reduces depth of getPrimeLoop recursion to logarithmic, so g++ can complete with default depth (without measurable time penalty):

constexpr int getPrimeLoop(int i, int k) {
    return (i == 0)?k:
        (i % 2)?getPrimeLoop(i-1, nextPrime(k + 1)):
        getPrimeLoop(i/2, getPrimeLoop(i/2, k));
}

Answer 2

I've given this a try myself, and written the following implementation:

template<unsigned... args> constexpr unsigned countArgs();
template<> constexpr unsigned countArgs() { return 0; }
template<unsigned head, unsigned... tail>
constexpr unsigned countArgs() { return 1 + countArgs<tail...>(); }

template<unsigned last>
constexpr unsigned lastArg() { return last; }
template<unsigned head, unsigned next, unsigned... tail>
constexpr unsigned lastArg() { return lastArg<next, tail...>(); }

template<unsigned i> constexpr bool isPrime() { return true; }
template<unsigned i, unsigned head, unsigned... tail>
constexpr bool isPrime()
{ return i % head && isPrime<i, tail...>(); }

template<bool found, unsigned i, unsigned... primesSoFar> struct nextPrime
{ static constexpr unsigned val =
    nextPrime<isPrime<i + 2, primesSoFar...>(), i + 2, primesSoFar...>::val; };
template<unsigned i, unsigned... primesSoFar> struct
nextPrime<true, i, primesSoFar...> { static constexpr unsigned val = i; };

template<unsigned n, unsigned... primesSoFar> struct nthPrimeImpl
{ static constexpr unsigned val = nthPrimeImpl<n - 1, primesSoFar...,
    nextPrime<false, lastArg<primesSoFar...>(), primesSoFar...>::val>::val; };
template<unsigned... primesSoFar> struct nthPrimeImpl<0, primesSoFar...>
{ static constexpr unsigned val = lastArg<primesSoFar...>(); };

template<unsigned n>
constexpr unsigned nthPrime() {
  return n == 1 ? 2 : nthPrimeImpl<n - 2, 3>::val;
}

constexpr unsigned p512 = nthPrime<512>();
static_assert(p512 == 3671, "computed incorrectly");

This requires increasing the maximum template depth of gcc to more than the default of 900 (in my gcc 4.7.2), e.g. by passing -ftemplate-depth=1200. And it is way too slow: it takes about 3 minutes on my hardware. So I very much hope for some more efficient code in a different answer.

In terms of computation method, the above does something like trial division. A sieve of Eratosthenes might perform better, but so far I couldn't think of a way to write that in a constexpr-compatible fashion.

Answer 3

The answer by Mike Kinghan got me thinking along lines I hadn't though before. If template instantiation is such a problem which causes such severe memory consumption, how can we reduce that? I eventually came up with a scheme where instead of an argument pack for all primes found so far, I use a chain of types, each one referencing the one before, and a chain of static functions in these types which can use the ones from types before.

The result I'll paste below is still a lot slower than the one suggested by zch, but I find it interesting enough to share, since it might be a useful approach for other applications.

template<unsigned N> struct NthPrime {
  typedef NthPrime<N - 1> previous;
  static constexpr unsigned prime = previous::nextPrime();
  static constexpr unsigned nextPrime() { return nextCoprime(prime + 2); }
  static constexpr unsigned nextCoprime(unsigned x) {
    // x is a candidate. We recurse to obtain a number which is
    // coprime to all smaller primes, then check that value against
    // the current prime.
    return checkCoprime(previous::nextCoprime(x));
  }
  static constexpr unsigned checkCoprime(unsigned x) {
    // if x is coprime to the current prime as well, then it is the
    // next prime. Otherwise we have to try the next candidate.
    return (x % prime) ? x : nextCoprime(x + 2);
  }
};

template<> struct NthPrime<1> {
  static constexpr unsigned prime = 2;
  static constexpr unsigned nextPrime() {
    return 3;
  }
  static constexpr unsigned nextCoprime(unsigned x) {
    return x; // x is guaranteed to be odd, so no need to check anything.
  }
};

template<unsigned n>
constexpr unsigned nthPrime() {
  return NthPrime<n>::prime;
}

constexpr unsigned p512 = nthPrime<512>();
static_assert(p512 == 3671, "computed incorrectly");

The beast above requires modifications both to the constexpr depth and the template depth. The following values are sharp bounds for my compiler.

time g++-4.7.2 -c -fconstexpr-depth=519 -ftemplate-depth=2042 -std=c++11 foo.cc

real    0m0.397s
user    0m0.368s
sys     0m0.025s

Answer 4

I doubt that your 1 second goal is within reach of any hardware that doesn't have security guards. However I believe the following meta-program can close down on it a lot:

#include <type_traits>

template<unsigned N>
using boolean = std::integral_constant<bool,N>;

template<unsigned N>
constexpr bool is_co_prime()
{
    return true;
};

template<unsigned N, unsigned D>
constexpr bool is_co_prime()
{
    return N % D != 0;
};

template<unsigned N, unsigned D0, unsigned D1, unsigned ...Di>
constexpr bool is_co_prime()
{
    typedef typename std::conditional<
        is_co_prime<N,D1>(),
        typename std::conditional<
            is_co_prime<N,Di...>(),
            boolean<is_co_prime<N,D0>()>,
            std::false_type
        >::type,
        std::false_type
    >::type type;
    return type::value;
};

template<unsigned N>
constexpr unsigned inc()
{
    return N == 2 ? 3 : N + 2;
}

template<unsigned Counter, unsigned Candidate, unsigned ...Primes>
struct nth_prime;

template<unsigned Candidate, unsigned Prime, unsigned ...Primes>
struct nth_prime<0,Candidate,Prime,Primes...>
{
    static const unsigned value = Prime;
};

template<unsigned Counter, unsigned Candidate = 2, unsigned ...Primes>
struct nth_prime
{
    typedef typename std::conditional<
        is_co_prime<Candidate,Primes...>(),
        nth_prime<Counter - 1,inc<Candidate>(),Candidate,Primes...>,
        nth_prime<Counter,inc<Candidate>(),Primes...>
    >::type type;
    static const unsigned value = type::value;
};

#include <iostream>

using namespace std;

int main()
{
    cout << nth_prime<512>::value << endl;
    return 0;
}

I'll call this meta-program MyNthPrime and call your one YourNthPrime.

Your seem to have much heftier hardware than mine and certainly more memory. I have a routine Lenovo ThinkPad T420, 4-core i5, 8GB RAM, 8GB swap, running Linux Mint 14, kernel 3.5.0. You report it takes you 3 mins. to build your YourNthPrime. Measured by the time command, it takes me 35 mins. to build YourNthPrime, with no apps running but the terminal and system monitor.

The compiler is GCC 4.7.2 and the commandline the options are:

-Wall -O0 -std=c++11 -ftemplate-depth=1200

That elapsed time breaks down as:

real    34m2.534s
user    3m40.414s
sys     0m33.450s

It takes me 1.5 mins to build MyNthPrime, with:

-Wall -O0 -std=c++11 -ftemplate-depth=2100

and the elapsed time breaks down as:

real    1m27.832s
user    1m22.205s
sys     0m2.612s

That -ftemplate-depth=2100 is not a transposition typo. More of this shortly.

MyNthPrime isn't fairly and squarely 23 times faster than YourNthPrime. The broken-down timings suggest that MyNthPrime is actually around 2.75 times as fast as as YourNthPrime on user-time. But they also show that YourNthPrime's build really lost the farm on real time. What was it doing for the fruitless 9/10ths of its existence? Swapping, of course.

Both builds scoffed 95% of my 8GB system RAM within 45 secs., but MyNthPrime topped out around there and didn't swap. YourNthPrime went on eating swap space up to a peak 3.9GB, and long before that all my CPUs were dozing.

This point is notable when you take it with the fact that MyNthPrime needs getting on for twice as much -ftemplate-depth as YourNthPrime. Folk wisdom is that an extravagant -ftemplate-depth is the road to ruin for a meta-program build time, because it equates to extravagant memory consumption, which only has to slide into heavy swapping and you're watching paint dry. But the run-off of YourNthPrime and MyNthPrime does not bear this out - quite the opposite. The lesson I take is that the depth to which you must instantiate recursive templates is not always a good measure of the quantity of template instantiating you must do, and it is the quantity that matters for your memory resources.

Although they don't look superficially similar, MyNthPrime and YourNthPrime, both implement the trial-division algorithm for prime generation. MyNthPrime is that much faster solely because it is rather more shrewdly coded to be sparing to of recursive template instantiations, and of the memory they gobble up.

YourNthPrime fields 4 recursive templates for the computation, all consuming the same recursive variadic template argument list. MyNthPrime fields 2: it's just giving the compiler about half as many enormous instantiations to do.

YourNthPrime (as I read it) perceives the potential efficiency of doing its trial-divisions in ascending order of the primes in hand - because the chances of successful division increase toward the smaller primes; and once a prime in hand exceeds 1/2 of the candidate number being divided, the chances are 0. Hit the most likely divisors first and you optimize your prospect of a quick verdict and exit. But an obstacle to exploiting this efficiency is the fact that the primes in hand are represented by a variadic template argument list with the largest always at its head. To surmount this obstacle YourNthPrime deploys the recursive variadic template function lastArg<>(), effectively to reverse the order in which the primes in hand are consumed in the divisions.

lastArg<>() presents the primes in hand to the template function:

template<unsigned i, unsigned head, unsigned... tail>
constexpr bool isPrime() {
    return i % head && isPrime<i, tail...>();
}

in ascending order for recursive trial division of the next candidate i by the primes head, tail.... It's here I think you look to lastArg<>() to pay its way by ensuring head is always he next-best prospect for getting you out of the costly righthand side of that &&.

But to achieve this lastArg<>() itself recursively traverses the whole list of primes in hand at each invocation, before you even get the opportunity for a verdict on i. So it would be cheaper just to let isPrime<>() traverse the primes in hand as far as it needs to, testing i as it goes, dispense with lastArg<>() and save all the recursive instantiations thereof.

The job done by isPrime<>() in YourNthPrime - the recursive trial division - is done in MyNthPrime by:

template<unsigned N, unsigned D0, unsigned D1, unsigned ...Di>
constexpr bool is_co_prime()
{
    typedef typename std::conditional<
        is_co_prime<N,D1>(),
        typename std::conditional<
            is_co_prime<N,Di...>(),
            boolean<is_co_prime<N,D0>()>,
            std::false_type
        >::type,
        std::false_type
    >::type type;
    return type::value;
};

is_co_prime<>() takes 10 lines to do what is_prime<>() does in one, I could have done it in one line also:

return is_co_prime<N,D0>() && is_co_prime<N,D1,Di...>();

could be the body of the function. But the ugly trounces the pretty for efficiency here. Every time is_prime<>() has to recurse into the tail, that tail is only one prime shorter than it was before. Every time is_co_prime<>() has to do the same thing the tail is two primes shorter than it was before. Its tail recursions are shallower in the worst case than those of is_prime<>() and there can only half as many.

is_co_prime<>() divides the candidate number N by the righthand - the smallest and likeliest - of any pair of available divisors first, and returns no-prime on success; otherwise recurses to any divisors still to the right of that pair and continues trial division of the next-but-one until success or exhaustion. Only at exhaustion does it resort to trial division by the larger of the the original pair - the least likely divisor of any it has tried. And likewise within each of the intervening smaller recursions, the least likely divisor gets tried last.

Though this algorithm can be seen to get speedily to the smaller and likelier divisors of N, an itch will persist to get straight to the smallest of them first and try them in true descending likelihood, as per lastArg<>(). We have to quell this itch with the recognition that any way of "getting straight to the smallest", when it is at the tail of the list, will be a meta-function that must itself recurse over the whole list before we even get started with trial divisions. The implementation of is_co_prime<>() gets us down the list "two steps at a time" doing the trial division as it does so.

True, sometimes it will unluckily "step over" the largest divisor of N first time and then won't find it again until and unless it bottoms out and recurses backwards up the list. But once N is big enough to make it matter there will mostly be at least one smaller divisor further on to the right and it would be really unlucky to miss all of them. So the risk of skipping the largest on the way down is no big deal. Remember too there aren't going to be any divisors on the way down till we reach the N/2 point. That means the worst-case waste of recursions by skipping the one-and-only divisor of some N on the way down is confined to the tail of the list from that point.

You speculate that a Sieve of Erathosthenes meta-program might compile faster, and you are right. As a prime generator, Sieve has better theoretical complexity than Trial Division. An elegant template meta-program implementation by Peter Simons, is here, dating from 2006 or before. (And as Peter Wood commented, a non-terminating Sieve of Erathosthenes meta-program broke the news that the C++ template system is Turing-complete.) With C++11 facilities Simons' meta-program can be much abbreviated but I don't think made much more efficient. Just it stands, Simons Sieve can generate at compiletime all the primes up to the 512th in under 9 secs. on my ThinkPad. It needs -ftemplate-depth=4000 or thereabouts to do it, but only about 0.5GB of memory - an even more arresting counterexample to the folk wisdom that big template-depth == big memory usage.

So Erathosthenes' Sieve leaves Trial Division floundering in the dust. Sadly, for the purpose of your exercise, Sieve is no use. It's called the sieve because we have to input an integer upper limit U and it sieves out of the composite numbers between 2 and U, leaving the primes. Thus to apply it for finding exactly the Nth prime = Pn and not possibly any other, you must anyhow compute Pn in some other way, give Sieve an upper limit of Pn + 1 (or Pn + 2, for Pn > 2), then throw away all the Pi, 2 <= Pi < Pn, that it returns to you, and keep just the Pn that you had already computed. This is a no-op.

Several of the commenters make the point that the identity any Nth prime you could possibly generate by compiletime meta-programming will be known beforehand or calculable beforehand by much simpler and vastly faster means. I can't disagree, but I support your general point that with C++11 facilities, TMP takes a huge stride toward real-world utility that is well worth exploring - the more so as whatever takes a minute to compile right now will take a second within in a decade.

Meanwhile, without even leaving our uncannily sophisticated C++ compilers, with TMP problems like this we can still experience the nature of programming an early computer, with clock speed and memory in the K's, in a "language" modelled tightly - but within arcane constraints! - on classical recursive function theory. That's why you've really gotta love them!