使用 Horner 方法进行多项式求值的 C++ constexpr

Posted 2023-02-22

【中文标题】使用 Horner 方法进行多项式求值的 C++ constexpr

【英文标题】：C++ constexpr for polynomial evaluation with Horner's method

【发布时间】：2018-05-08 11:46:31

【问题描述】：

我希望能够使用霍纳方法评估多项式的​​导数，并将结果用作constexpr。这看起来非常平凡，但我错过了一些明显的东西，因为编译器说我超过了最大递归深度。核心递归发生在这里：

template<size_t d, size_t i, typename C, typename X>
constexpr X evalImpl(const C &c, const X &x) 
    return i >= (C::SizeAtCompileTime - 1 - d) ? 1 : evalImpl<d, i + 1, C, X>(c, x);

你应该不需要知道Horner的方法就知道这里发生了什么，所以我已经尽可能地剥离了代码，甚至去掉了x的使用方式，因为它似乎与问题无关我有。

这个想法是，当索引i 等于多项式Degree<C>::value 的次数减去导数d 的阶数时，递归应该停止。否则，它应该增加索引i 并重试。

我正在使用表单调用来调用上述递归

 eval<derivative, 0>(c, x)

其中c 是Eigen::Matrix<double,1,7> 类型的特征矩阵，x 是双精度数。这个想法是从 0 开始，并从本质上算到多项式的次数。

编译器错误消息的格式为

In file included from /mnt/c/proj/src/main.cpp:11:0:
/mnt/c/proj/src/polynomial.h: In instantiation of 'constexpr X anonymous::evalImpl(const C&, const X&) [with long unsigned int d = 1ul; long unsigned int i = 898ul; C = Eigen::Matrix<double, 1, 7>; X = double]':
/mnt/c/proj/src/polynomial.h:74:108:   recursively required from 'constexpr X anonymous::evalImpl(const C&, const X&) [with long unsigned int d = 1ul; long unsigned int i = 1ul; C = Eigen::Matrix<double, 1, 7>; X = double]'
/mnt/c/proj/src/polynomial.h:74:108:   required from 'constexpr X anonymous::evalImpl(const C&, const X&) [with long unsigned int d = 1ul; long unsigned int i = 0ul; C = Eigen::Matrix<double, 1, 7>; X = double]'
/mnt/c/proj/src/polynomial.h:109:39:   required from 'constexpr X Polynomial::eval(const C&, const X&) [with long unsigned int d = 1ul; C = Eigen::Matrix<double, 1, 7>; X = double]'
/mnt/c/proj/src/main.cpp:306:66:   required from here
/mnt/c/proj/src/polynomial.h:74:108: fatal error: template instantiation depth exceeds maximum of 900 (use -ftemplate-depth= to increase the maximum)
         return i >= (C::SizeAtCompileTime - 1 - d) ? 1 : evalImpl<d, i + 1, C, X>(c, x);

【问题讨论】：

我不关注。多项式的次数为 7。这不应超过 7 次

一般来说，这段代码不会用于生成编译时表达式。大多数时候，我会将它与双打一起使用。但是，在某些情况下，我们需要将输出作为模板参数

【参考方案1】：

这里的条件：

return i >= (C::SizeAtCompileTime - 1 - d) ? 1 : evalImpl<d, i + 1, C, X>(c, x);

不是if constexpr。因此无论i >= (C::SizeAtCompileTime - 1 - d) 是true 还是false，剩下的总是会被实例化。因此递归不会随心所欲地停止。

改成这样：

if constexpr (i >= (C::SizeAtCompileTime - 1 - d)) 
    return 1;
 else  
    return evalImpl<d, i + 1, C, X>(c, x);  

编辑：

如果你无法访问 C++17，请使用标签调度：

template<size_t d, size_t i, typename C, typename X>
constexpr X evalImpl_impl(const C &c, const X &x, std::true_type) 
    return 1;


template<size_t d, size_t i, typename C, typename X>
constexpr X evalImpl_impl(const C &c, const X &x, std::false_type) 
    return evalImpl_impl<d, i + 1, C, X>(c, x, std::integral_constant<bool, C::SizeAtCompileTime-1-d <= i+1>);


template<size_t d, size_t i, typename C, typename X>
constexpr X evalImpl(const C &c, const X &x) 
    return evalImpl_impl(c, x, std::integral_constant<bool, C::SizeAtCompileTime-1-d <= i>);

【讨论】：

三元运算符不是 constexpr？

@bremen_matt 2 个分支都应该被实例化。短路只影响运行时评估，不影响实例化。

这是 C++11 吗？还是 14 岁？

@bremen_matt if constexpr 是 C++17。您使用三元运算符的方法不适用于任何版本。

我一直在尝试使用模板专业化将其分解为两个函数，一个是停止条件，但没有成功。我想我将不得不切换回这种方法。我需要 C++11