在编译时计算第n个素数

Question

constexpr在我看来,C++ 11的功能,以及模板参数包应该足够强大,可以执行一些相当复杂的计算.我有一个实际应用的一个可能的例子是在编译时计算第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深度,但它很容易适应时间:

$ 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内存,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的构建确实在实时失去了农场.它为存在的十分之九的无果而做了什么？当然,交换.

两个版本都在45秒内嘲笑了我的8GB系统RAM的95%,但MyNthPrime 在那里突然出现并没有交换.YourNthPrime继续吃交换空间达到3.9GB的峰值,很久以前我的所有CPU都在打瞌睡.

当你把它的事实,这一点是值得注意的MyNthPrime 需要获得关于两倍-ftemplate-depth的YourNthPrime.民间的智慧是,奢侈-ftemplate-depth是毁灭元程序构建时间的道路,因为它等同于奢侈的记忆消耗,只需要进行大量交换,你就会看到油漆干燥.但是YouNthPrime和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由右手-最小的,有可能的-任何一对可用的除数的第一,并成功返回不贷; 否则,对任何仍然在该对权利之外的任何除数进行递归,并继续对下一个但是继续进行试验,直到成功或用尽为止.只有在筋疲力尽的情况下才能通过原始对中较大的一对来进行试验分裂 - 这是它所尝试过的最不可能的除数.同样在每个中间较小的递归中,最不可能的除数最后被尝试.

虽然可以看出这个算法可以快速地获得更小和更有可能的除数N,但是一个痒将继续直接到达它们中最小的一个并按照真正的下降可能性来尝试它们lastArg<>().我们不得不平息这种痒,认识到任何"直到最小"的方式,当它在列表的尾部时,将是一个元函数,在我们开始之前必须自己递归整个列表有审判部门.实施is_co_prime<>()让我们在列表中"一步两步"进行试验分工.

没错,有时它会不幸地"跨过" N 第一次最大的除数,然后再也找不到它,直到它除非最低点并向后递回列表.但是一旦N大到足以使其变得重要,那么在右边进一步将至少有一个较小的除数,错过所有这些将是非常不幸的.因此,在下降途中跳过最大的风险并不是什么大问题.记住 ,在我们到达目的地之前,不会有任何除数N/2.这意味着通过N在向下的过程中跳过某些的一个和唯一的除数而最坏情况下的递归浪费被限制在该点的列表的尾部.

你推测,Erathosthenes元程序的Sieve可能编译速度更快,你是对的.作为主要发电机,Sieve具有比试验部更好的理论复杂性.通过优雅的模板元编程实现彼得·西蒙斯,就是 在这里,从2006年或之前约会.(正如彼得伍德评论的那样,一个非终结的Erathosthenes元程序筛选器爆出了C++模板系统图灵完备的消息.)使用C++ 11设备,Simons的元程序可以缩写很多,但我不知道认为效率更高.就这一点而言,Simons Sieve可以在编译时生成所有素数,直到第512个不到9秒.在我的ThinkPad上.它需要-ftemplate-depth=4000或者需要这样做,但只有大约0.5GB的内存 - 这是民间智慧的一个更具吸引力的反例,即大模板深度==大内存使用率.

所以Erathosthenes的Sieve让审判部在尘土中挣扎.可悲的是,为了锻炼,Sieve是没用的.它被称为筛子,因为我们必须输入一个整数上限,U并从2和2之间的复合数中筛选U,留下素数.因此,要应用它来准确找到Nth prime = Pn而不是任何其他的,你必须Pn以其他方式计算,给Sieve一个上限Pn + 1(或Pn + 2,for Pn > 2),然后扔掉Pi, 2 <= Pi < Pn它返回给你的所有,并保持Pn你已经计算过的.这是一个无操作.

一些评论者指出,您可能通过编译时元编程生成的任何Nth prime的身份将事先知道或通过更简单和更快速的方式预先计算.我不能不同意,但我支持你的一般观点,即使用C++ 11工具,TMP向现实世界的实用工具迈出了巨大的一步,这是非常值得探索的 - 现在需要花一分钟时间编译就更需要十年之内的第二次.

与此同时,即使离开了我们不可思议的复杂C++编译器,如果遇到这样的TMP问题,我们仍然可以体验编程早期计算机的本质,即K中的时钟速度和内存,紧密模拟的"语言" - 但是在神秘的约束下! - 关于经典递归函数理论.这就是为什么你真的要爱他们!