我正在尝试实现高斯-勒让德求积，我想要一个模板化的函数将点数作为模板参数。现在我有这个：
template<int number_of_quadrature_points>
double gaussian_quadrature_integral_core(double (*f)(double), double from, double to){
double scaling_factor = (to-from)/2;
double average_factor = (from+to)/2;
std::array<double, number_of_quadrature_points> factors;
std::array<double, number_of_quadrature_points> points;
if constexpr(number_of_quadrature_points == 2){
factors = {1, 1};
points = {-1.0/sqrt(3), 1.0/sqrt(3)};
}
if constexpr(number_of_quadrature_points == 3){
factors = {5.0/9.0, 8.0/9.0, 5.0/9.0};
points = {-sqrt(3.0/5.0), 0, sqrt(3.0/5.0)};
}
double sum = 0;
for(int i = 0; i < number_of_quadrature_points; i++){
sum += factors.at(i)*((*f)(scaling_factor*points.at(i)+average_factor));
}
sum *= scaling_factor;
return sum;
}
正如您所看到的，当模板参数发生变化时，不仅数组大小发生变化，内容也发生变化，但对于给定的大小，内容是众所周知的。出于这个原因，我认为如果std:：数组是const-static会更好，因为函数被调用了很多次。
目前，我只使用了if constexpr来声明数组，但我如何使用它来定义和声明数组，使其在if constexpr范围之外可见，并且数组只定义了一次？
添加两个辅助函数就足够了(如果您使用C++20(：
template<unsigned N>
constexpr auto init_factors() {
std::array<double, N> rv;
if constexpr(N == 2){
rv = {1., 1.};
} else {
rv = {5.0/9.0, 8.0/9.0, 5.0/9.0};
}
return rv;
}
template<unsigned N>
constexpr auto init_points() {
std::array<double, N> rv;
if constexpr(N == 2){
rv = {-1.0/std::sqrt(3.), 1.0/std::sqrt(3.)};
} else {
rv = {-std::sqrt(3.0/5.0), 0, std::sqrt(3.0/5.0)};
}
return rv;
}
template<unsigned number_of_quadrature_points>
double gaussian_quadrature_integral_core(double (*f)(double), double from,
double to)
{
static constexpr auto factors = init_factors<number_of_quadrature_points>();
static constexpr auto points = init_points<number_of_quadrature_points>();
[...]
为了防止使用错误的点数，您可以添加一个static_assert
template<unsigned number_of_quadrature_points>
double
gaussian_quadrature_integral_core(double (*f)(double), double from,
double to)
{
static_assert(number_of_quadrature_points==2||number_of_quadrature_points==3);
或者如果你想稍后进行专业化，请使用SFINAE阻止匹配：
#include <type_traits>
template<unsigned number_of_quadrature_points>
std::enable_if_t<number_of_quadrature_points==2||number_of_quadrature_points==3,
double>
gaussian_quadrature_integral_core(double (*f)(double), double from,
double to)
{
您可能有模板变量：
template <std::size_t N>
static constexpr std::array<double, N> factors;
template <std::size_t N>
static constexpr std::array<double, N> points;
template <>
constexpr std::array<double, 2> factors<2>{{1, 1}};
template <>
constexpr std::array<double, 2> points<2>{{-1.0 / sqrt(3), 1.0 / sqrt(3)}};
template <>
constexpr std::array<double, 3> factors<3>{{5.0 / 9.0, 8.0 / 9.0, 5.0 / 9.0}};
template <>
constexpr std::array<double, 3> points<3>{{-sqrt(3.0 / 5.0), 0, sqrt(3.0 / 5.0)}};
然后
template<int number_of_quadrature_points>
double gaussian_quadrature_integral_core(double (*f)(double), double from, double to)
{
const double scaling_factor = (to - from) / 2;
const double average_factor = (from + to) / 2;
double sum = 0;
for(int i = 0; i < number_of_quadrature_points; i++){
sum += factors<number_of_quadrature_points>[i]
* ((*f)(scaling_factor * points<number_of_quadrature_points>[i] + average_factor));
}
sum *= scaling_factor;
return sum;
}
演示
请注意，如果没有constexprsqrt(std::不是(，则必须将constexpr替换为const。
您可以在本主题中使用类似的内容：有没有一种方法可以在C++模板专用化中对常值参数设置条件？
因此，我们使用std:：enable_if和SFINAE创建了两个模板专业化。我们通过模板参数number_of_quadrature_points来区分它们。这样，我们就有了全局参数，而不必多次定义和实例化。此代码使用c++17进行编译。
此外，我还建议使用std:：function<gt；而不是指向函数的指针。
#include <array>
#include <cmath>
#include <iostream>
#include <functional>
template<int number_of_quadrature_points, typename E=void>
struct gaussian_quadrature_params
{
};
template<int number_of_quadrature_points>
struct gaussian_quadrature_params<number_of_quadrature_points, std::enable_if_t<(number_of_quadrature_points==2)> >
{
constexpr static const std::array<double, number_of_quadrature_points> factors = {1, 1};
constexpr static const std::array<double, number_of_quadrature_points> points = {-1.0/sqrt(3), 1.0/sqrt(3)};
};
template<int number_of_quadrature_points>
struct gaussian_quadrature_params<number_of_quadrature_points, std::enable_if_t<(number_of_quadrature_points==3)> >
{
constexpr static const std::array<double, number_of_quadrature_points> factors = {5.0/9.0, 8.0/9.0, 5.0/9.0};
constexpr static const std::array<double, number_of_quadrature_points> points = {-sqrt(3.0/5.0), 0, sqrt(3.0/5.0)};
};
double f(double x)
{
return x;
}
template<int number_of_quadrature_points>
double gaussian_quadrature_integral_core(std::function<double(double)> f, double from, double to){
double scaling_factor = (to-from)/2;
double average_factor = (from+to)/2;
double sum = 0;
for(int i = 0; i < number_of_quadrature_points; i++){
sum += gaussian_quadrature_params<number_of_quadrature_points>::factors.at(i)*(f(scaling_factor*gaussian_quadrature_params<number_of_quadrature_points>::points.at(i)+average_factor));
}
sum *= scaling_factor;
return sum;
}
int main()
{
std::cout << gaussian_quadrature_integral_core<2>(f, -1.0, 1.0) << std::endl;
std::cout << gaussian_quadrature_integral_core<3>(f, -1.0, 1.0) << std::endl;
}