Define and declare a const std::array with constexpr if

Question

I am trying to implement Gauss-Legendre quadrature, and I want a templated function that takes the number of points as a template parameter. Right now I have this:

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

As you see, when the template parameter changes, not only the array size changes, but also the content is changed, but for a given size the content is well known. For this reason I think it would be better if the std::arrays were const static, as the function is called many times.

Right now I have only managed to use if constexpr to declare the array, but how can I use it both to define and declare the array so It is visible outside the if constexpr scope and the arrays are only defined once?

What happens if someone calls your function with N different from 2 or 3? if constexpr are great, but I would go for a specialized template in that case (for the factors/points selection). — Cedric
– Cedric, Commented Dec 17, 2020 at 14:04
@Cedric this is indeed a problem, right now I just wanted to do a static assert, but I will reconsider it as Your comment made me aware such a solution might not be elegant. — Karol Szustakowski
– Karol Szustakowski, Commented Dec 17, 2020 at 14:08

Ted Lyngmo · Accepted Answer · 2020-12-17 14:13:01Z

1

Adding two helper functions could be enough (if you're using 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>();
[...]

To prevent usage with the wrong number of points, you could add a 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);

...or prevent matching using SFINAE if you want to make a specialization later:

#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)
{

edited Dec 17, 2020 at 14:13

answered Dec 17, 2020 at 13:57

Ted Lyngmo

124k7 gold badges93 silver badges155 bronze badges

Sign up to request clarification or add additional context in comments.

2 Comments

Karol Szustakowski Over a year ago

Something like this also came to my mind, but I was thinking if it was possible without the helper functions. It is still a good solution though, thank you.

Passer By Over a year ago

@KarolSzustakowski Use lambdas if you prefer them.

Jarod42 · Accepted Answer · 2020-12-17 15:46:26Z

1

You might have template variable:

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

and then

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

Demo

Notice that you have to replace constexpr by const if you don't have constexpr sqrt (which std:: isn't).

edited Dec 17, 2020 at 15:46

answered Dec 17, 2020 at 15:32

Jarod42

227k15 gold badges214 silver badges350 bronze badges

1 Comment

Ted Lyngmo Over a year ago

... and it works even in C++14, which is nice.

trebuhg · Accepted Answer · 2020-12-17 14:38:35Z

You can use something like in this topic: is there a way to put condition on constant value parameter in C++ template specialization?

So we are creating two template specialisations using std::enable_if and SFINAE. We differentiate them by template parameter number_of_quadrature_points. This way we have global parameters which does not have to be defined and instantiated multiple times. This code compiles using c++17.

Also I suggest to use modern approach with std::function<> instead of pointer to function.

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

You might simply have specialization for template<> struct gaussian_quadrature_params<2> with enabler (you can keep enabler to allow real formula (as odd numbers greater than 42)).

Ping Chia · Accepted Answer · 2023-12-21 01:44:47Z

How about

// N: number_of_quadrature_points
template<int N>
double gaussian_quadrature_integral_core(double (*f)(double), double from, double to)
{
    constexpr std::array<double, N> factors = []()
        ->std::array<double, N>{
        if constexpr(N == 2)
            return {1.0, 1.0};
        else if constexpr(N == 3)
            return {5.0 / 9.0, 8.0 / 9.0, 5.0 / 9.0};
        // ... other N cases
    }();

    constexpr std::array<double, N> points= []()->auto{
        if constexpr(N == 2)
            return std::array<double, N>{-1.0 / std::sqrt(3), 1.0 / std::sqrt(3)};
        else if constexpr(N == 3)
            return std::array<double, N>{-std::sqrt(3.0 / 5.0), 0, std::sqrt(3.0 / 5.0)};
        // ... other N cases
    }();

    double scaling_factor = (to - from) / 2;
    double average_factor = (from + to) / 2;
    double sum = 0;
    for (int i = 0; i < N; i++)
        sum += factors.at(i)*((*f)(scaling_factor * points.at(i) + average_factor));

    sum *= scaling_factor;
    return sum;
}

which declare and define the arrays using if constexpr.

Collectives™ on Stack Overflow

Define and declare a const std::array with constexpr if

4 Answers 4

2 Comments

1 Comment

1 Comment

Comments

Your Answer

Linked

Hot Network Questions

Collectives™ on Stack Overflow

4 Answers 4

2 Comments

1 Comment

1 Comment

Comments

Your Answer

Sign up or log in

Post as a guest

Linked

Related