17

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

我写的代码做得很好:

#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 9错误而失败

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 和深度限制为 512 的编译字符串constexpr。(我的 g++ 的默认值)。它应该适用于所有正 int32 或至少适用于高达 10 9 +9的数字

4

6 回答 6

12

您可以使用分治算法将范围分成两半,从而将空间要求从线性更改为对数。此方法使用分而治之,并且仅测试奇数因素(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 的答案中窃取缩进样式。

于 2013-08-19T05:05:04.283 回答
7

这是我的尝试。使用Miller-Rabin 素数测试constexpr的确定性变体,用于最多 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;
}
于 2013-08-19T06:11:38.013 回答
4

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

更新:修复了 Newton-Raphson 整数平方根

这段代码不是最理想的——将所有测试除以偶数(甚至可能是三的倍数)显然会加快编译时间——但它可以工作,即使使用大约 10 10 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 { };
于 2013-08-19T07:24:14.607 回答
3

你可以看看 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
于 2013-08-18T21:10:29.020 回答
2

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

我可以想到两种增加上限的策略:

  1. 改进算法。如果只检查奇数因子,则将迭代次数减半。上限上升了四倍。这仍然没有接近 10 亿,但您当然可以通过更接近理想筛子来实现。

  2. 使用typedef声明来预先评估部分元函数,并依靠编译器的记忆策略来阻止该部分被全面重新评估。

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

于 2013-08-19T02:07:37.630 回答
-1

没有 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 };
};
于 2015-11-25T09:04:38.223 回答