编译时素数检查答案

【问题标题】：Compile time prime checking编译时素数检查
【发布时间】：2013-08-20 15:35:31
【问题描述】：

我需要检查编译时是否有一些整数素数（将布尔值作为模板参数）。

我写的代码做得很好：

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

它适用于 ~1000000 的数字，但失败并出现错误 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

我无法增加深度限制。有没有可能减少我使用的深度？

我想要实现的目标：我需要检查编译时是否为常数素数无需更改编译字符串，模板深度限制为 900，constexpr 深度限制为 512。（我的 g++ 的默认设置）。它应该适用于所有正 int32 或至少适用于高达 10⁹+9

的数字

【问题讨论】：

为什么不使用 constexpr？看看这里：cpptruths.blogspot.no/2011/07/…
所以，你被 C++ 大师控制了，带着... 作业... :)
@olevegard 我不知道，但并非所有声称支持 C++11 的编译器都有 constexpr。（我在看你VS2012...）
@olevegard，嘿，我只是忘记了它的存在。您可以将其发布为答案吗？
@sehe，信不信由你，我确实自己编写了这段代码，这是我自己的想法。我被我控制了吗？可能是。但我不会说他（我）是大师

标签： c++ templates c++11 template-meta-programming compile-time

【解决方案1】：

您可以通过使用分治算法将范围分成两半，将空间要求从线性更改为对数。该方法采用分治法，只测试奇数因子（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))));
}

编辑：使用编译时 sqrt 将搜索空间绑定到上限 (sqrt(n))，而不是 n / 2。现在可以根据需要在 Coliru 上计算 is_prime(100000007)（以及 is_prime(1000000000039ULL)）而不会爆炸.

对于糟糕的格式表示歉意，我仍然没有为 C++11 折磨的 constexpr 子语言找到一种舒适的风格。

编辑：清理代码：用另一个函数替换宏，将实现细节移动到细节命名空间中，从 Pablo 的答案中窃取缩进样式。

【讨论】：

“为可怕的格式道歉......” 别担心，它比等效的模板元编程格式好数千倍:)
我喜欢这样。我也用它来检查偶数，它看起来更干净一些。谢谢。

【解决方案2】：

这是我的尝试。将 constexpr 和 deterministic variant of the Miller-Rabin primality test 用于最多 4,759,123,141 的数字（应该涵盖所有 uint32，但您可以轻松更改底漆检查器集以覆盖更大的范围）

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

【讨论】：

【解决方案3】：

constexpr 可能更容易处理，但使用纯模板实例化并没有真正的问题。

更新：修正 Newton-Raphson 整数平方根

这段代码不是最理想的——将所有测试除以偶数（甚至可能是三的倍数）显然会加快编译时间——但它可以工作，即使素数约为 10^{10 sup> gcc 使用不到 1GB 的 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 { };

【讨论】：

【解决方案4】：

你可以看看 constexpr。它的语法比模板元编程更友好（至少如果你不熟悉像我这样的模板。）你不能使用 if 或任何循环。但是使用递归和三元运算符，您可以使用模板元编程做几乎所有可以做的事情，而且它通常运行得更快。

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

这是一个使用在线编译器的工作示例： http://coliru.stacked-crooked.com/view?id=6bc10e71b8606dd2980c0c5dd982a3c0-6fbdb8a7476ab90c2bd2503cd4005881

由于它是在编译时执行的，你可以做一个静态断言来测试它。

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

如果x 不是素数，则断言将失败，这意味着编译将失败。如果 x 是素数，编译将成功，但不会生成任何输出。

如果你想检查非常大的数字，你可以增加 constexpr 深度

-fconstexpr-depth=930000

我还没有测试过它支持多大的数字，但我认为它因编译器而异。

如果你想自己测试一下：

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

编译

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

【讨论】：

嗯，我注意到它的深度限制在（甚至限制更低）：ideone.com/VGJxnr
阅读文章，您可以使用-fconstexpr-depth=9300 设置限制如果您手动设置深度，constexpr 可以比模板元编程“更深入”。至少根据文章，我自己没试过。
我无法更改编译字符串。如果可以的话，我会为自己的代码更改它。
coliru.stacked-crooked.com 它为我使用这个编译器工作：coliru.stacked-crooked.com/…
之所以有效，是因为 7 是一个非常小的数字。例如，它不会为 1000000007 编译

【解决方案5】：

从语言的角度来看，解决方案是增加深度限制。该程序是正确的，只是它需要“太多”的迭代。但是你已经说过你不想增加它。（看起来所需的模板深度是 (sqrt(N) + C)，其中 C 是一个小常数，因此对于系统上的最大深度 900，您当前的实现将工作到 810000。）

我可以想到两种策略来增加上限：

改进算法。如果只检查奇数因子，则将迭代次数减半。上限上升了四倍。这仍然没有接近 10 亿，但您当然可以通过更接近理想筛子来达到目标。
使用typedef 声明来预先评估部分元函数，并依靠编译器的记忆策略来阻止该部分被全面重新评估。

此策略不适用于严重依赖先前迭代结果的元函数，但在您的情况下，您可以检查最后 900 个因子，然后检查最后 1800 个因子将自动使用结果的缓存副本从最后的 900 开始。这不是 C++ 标准规定的，严格来说是不可移植的，但另一方面也没有关于这些递归限制的任何内容。

【讨论】：

【解决方案6】：

没有 constexpr 的 C++，IsPrime::Value 给出编译时结果。诀窍是反复尝试除以 i=3,5,7,... 直到 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 };
};

【讨论】：

嗯，它与OP中的基本相同，并没有解决想要的问题ideone.com/u45RNt
枚举总是在编译时进行评估，例如 IsPrime::Value。你在寻找不同的东西吗？
我需要它来处理大数字而不达到实例化深度的限制。