c++ - 在编译时计算第 n 个素数

Question

在我看来，C++11 特性constexpr以及模板参数包应该足够强大，可以执行一些相当复杂的计算。我有一个实际应用的一个可能示例是在编译时计算第 n 个素数。

我正在寻求实现这种计算的方法。如果提出了多个解决方案，比较它们可能会很有趣。

为了让您了解我的性能期望：我希望在合理的桌面硬件上，一些代码可以在不到一秒的编译时间内找到第 512 个素数（即 3671）。

Answer 1

我实现了最简单的方法，根本不使用模板，它可以工作：

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");

它需要增加 constexpr-depth 一点，但它很容易适应时间：

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

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

以下技巧将getPrimeLoop递归深度减少到对数，因此 g++ 可以使用默认深度完成（没有可测量的时间损失）：

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

我怀疑您的 1 秒目标是否在任何没有安全防护的硬件的范围内。但是我相信以下元程序可以关闭它很多：

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

我将把这个元程序称为 MyNthPrime并将你的一个称为 YourNthPrime。

你的硬件似乎比我的要重得多，当然内存也更多。我有一个常规的联想 ThinkPad T420，4 核 i5，8GB RAM，8GB 交换，运行 Linux Mint 14，内核 3.5.0。你报告它需要你 3 分钟。构建您的YourNthPrime。按time命令测量，我需要 35 分钟。构建YourNthPrime，除了终端和系统监视器之外没有运行任何应用程序。

编译器是 GCC 4.7.2，命令行选项是：

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

经过的时间分解为：

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

构建MyNthPrime需要 1.5 分钟，其中：

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

并且经过的时间分解为：

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

这-ftemplate-depth=2100不是换位错字。很快就会有更多。

MyNthPrime并不比YourNthPrime快 23 倍。分解的时间表明MyNthPrime在用户时间上实际上是YourNthPrime的 2.75 倍左右。但它们也表明YourNthPrime的构建确实实时丢失了农场。它为它存在的 9/10 的无结果做了什么？交换，当然。

两种构建都在 45 秒内嘲笑了我 8GB 系统 RAM 的 95%，但MyNthPrime 在那里达到顶峰并且没有交换。YourNthPrime继续占用高达 3.9GB 的交换空间，而在那之前，我所有的 CPU 都在打瞌睡。

当您考虑到MyNthPrime 需要两倍于-ftemplate-depthYourNthPrime的事实时，这一点是值得注意的。民间智慧是，奢侈-ftemplate-depth是元程序构建时间的毁灭之路，因为它等同于奢侈的内存消耗，它只需要滑入大量交换并且你正在看着油漆干涸。但是YourNthPrime和MyNthPrime的决胜并没有证明这一点 - 恰恰相反。我得到的教训是，你必须实例化递归模板的深度并不总是一个很好的数量衡量标准你必须做的模板实例化，它的数量对你的内存资源很重要。

虽然它们看起来并不相似，但 MyNthPrime和YourNthPrime都实现了质数生成的试除算法。MyNthPrime 的速度要快得多，仅仅是因为它的编码更加精明，可以节省递归模板实例化以及它们吞噬的内存。

YourNthPrime为计算提供 4 个递归模板，都使用相同的递归可变参数模板参数列表。 MyNthPrime字段 2：它只是为编译器提供了大约一半的巨大实例化。

YourNthPrime（正如我所读的那样）感知到按手头素数升序进行试除的潜在效率 - 因为成功除法的机会向较小的素数增加；一旦手中的素数超过被除数的 1/2，机会为 0。首先击中最有可能的除数，然后优化快速判断和退出的前景。但是利用这种效率的一个障碍是，手头的素数由可变参数模板参数列表表示，最大的 总是在它的头部。为了克服这个障碍， YourNthPrime部署了递归可变参数模板函数lastArg<>()，有效地反转了手头素数在除法中的消耗顺序。

lastArg<>()将手中的素数呈现给模板函数：

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

按升序对下一个候选i 除以素数进行递归试除head, tail...。在这里，我认为你希望lastArg<>() 通过确保head他永远是让你摆脱代价高昂的右手边的下一个最佳前景来付出代价&&。

但要实现这一点，lastArg<>()它本身会在每次调用时递归遍历手头的整个素数列表，甚至在您有机会对i. isPrime<>() 因此，只需要尽可能地遍历手头的素数，进行测试i，省去lastArg<>()并保存其所有递归实例，就会更便宜。

isPrime<>()在YourNthPrime中完成的工作- 递归试验部门 - 在MyNthPrime中由以下人员完成：

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<>()需要 10 行来完成is_prime<>()一项工作，我也可以在一行中完成：

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

可能是函数的主体。但是在这里，丑陋的东西在效率上胜过漂亮的东西。每次is_prime<>()都必须递归到尾部，该尾部仅比以前短一个素数。每次 is_co_prime<>()都必须做同样的事情，尾巴比以前短了两个素数。它的尾递归在最坏的情况下比那些更浅，is_prime<>()并且只有一半。

is_co_prime<>()将候选数除以N任何一对可用除数的右手边（最小和最可能的）first，并在成功时返回无素数；否则递归到仍然在该对右侧的任何除数，并继续对下一个除数进行尝试除法，直到成功或用尽。只有在用尽时，它才会诉诸原始对中较大的除数 - 它尝试过的任何除数中最不可能的除数。同样，在每个中间较小的递归中，最不可能的除数最后被尝试。

尽管可以看到该算法可以快速找到 的较小和更可能的除数N，但仍然会坚持首先直接找到其中最小的除数，并按照真实的降序尝试它们lastArg<>()。我们必须通过认识到任何“直接到最小”的方式，当它位于列表的尾部时，将是一个元函数，它必须在我们开始之前在整个列表上递归。与审判部门。的实施is_co_prime<>()使我们在“一次两步”的情况下进行试验划分。

诚然，有时它会不幸地“越过”N 第一次的最大除数，然后不会再次找到它，除非它触底并在列表中向后递归。但是一旦N大到足以让它变得重要，在右边大部分将至少有一个更小的除数，错过所有这些真的很不幸。因此，在下降过程中跳过最大的风险没什么大不了的。 还要记住，在我们到达终点之前，在下降的过程中不会有任何除数。N/2这意味着通过在向下的过程中跳过一些唯一除数来浪费递归的最坏情况N仅限于从该点开始的列表尾部。

您推测 Sieve of Erathothenes 元程序可能编译得更快，您是对的。作为素数生成器，Sieve 比 Trial Division 具有更好的理论复杂性。Peter Simons的 一个优雅的模板元程序实现在这里，可以追溯到 2006 年或更早。（正如 Peter Wood 评论的那样，一个非终止的 Erathothenes 元程序筛子爆料说 C++ 模板系统是图灵完备的。）有了 C++11 设施，Simons 的元程序可以大大简化，但我不知道。 t认为变得更有效率。就目前而言，Simons Sieve 可以在编译时在 9 秒内生成直到第 512 个素数的所有素数。在我的 ThinkPad 上。它需要-ftemplate-depth=4000或者差不多要这样做，但只有大约 0.5GB 的内存 - 一个更引人注目的反例，与大模板深度 == 大内存使用的民间智慧相比。

因此，埃拉托色尼的筛子让审判部门在尘埃中挣扎。可悲的是，为了您的锻炼目的，Sieve 没有用。它被称为筛子，因为我们必须输入一个整数上限U，它会从 2 和 之间的合数中筛出U，留下素数。因此，要将其应用于准确地找到第 N 个素数 =Pn而不可能是其他任何素数，无论如何您必须Pn以其他方式计算，给 Sieve 一个上限Pn + 1（或Pn + 2, for Pn > 2），然后丢弃所有Pi, 2 <= Pi < Pn返回给您的 ,并只保留Pn你已经计算过的。这是一个无操作。

一些评论者指出，您可能通过编译时元编程生成的任何第 N 个素数的身份都将是预先知道的，或者可以通过更简单和更快的方式预先计算。我不能不同意，但我支持您的一般观点，即使用 C++11 工具，TMP 朝着非常值得探索的实际实用程序迈出了一大步 - 更何况现在需要一分钟的时间来编译十年内的第二次。

同时，即使不离开我们极其复杂的 C++ 编译器，面对这样的 TMP 问题，我们仍然可以体验到对早期计算机进行编程的本质，其时钟速度和内存在 K 中，采用紧密建模的“语言”——但在神秘的约束范围内！- 关于经典递归函数理论。这就是为什么你真的要爱他们！

Answer 3

我自己尝试过，并编写了以下实现：

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");

这需要将最大模板深度gcc增加到超过默认值 900（在我的 gcc 4.7.2 中），例如通过传递-ftemplate-depth=1200. 而且太慢了：在我的硬件上大约需要 3分钟。所以我非常希望在不同的答案中有一些更有效的代码。

在计算方法方面，上面做了类似试验除法的事情。Eratosthenes的筛子可能会表现得更好，但到目前为止，我想不出一种方法来以与 constexpr 兼容的方式编写它。

Answer 4

Mike Kinghan的回答让我按照我以前从未想过的思路思考。如果模板实例化是一个导致如此严重的内存消耗的问题，我们该如何减少呢？我最终想出了一个方案，在这个方案中，我使用了一系列类型，每个类型都引用了之前的类型，而不是迄今为止找到的所有素数的参数包，以及这些类型中的静态函数链，这些类型可以使用类型中的那些前。

我将在下面粘贴的结果仍然比zch 建议的要慢很多，但我发现它很有趣，可以分享，因为它可能是其他应用程序的有用方法。

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");

上面的野兽需要修改 constexpr 深度和模板深度。以下值是我的编译器的明确界限。

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