Compile-time Lagrange polynomials in C++

Question

I have implemented Lagrange polynomials, i.e.

\$l_j(x) := \prod_{\begin{smallmatrix}0\le m\le k\\ m\neq j\end{smallmatrix}}\frac{x-x_m}{x_j-x_m} = \frac{(x-x_0)}{(x_j-x_0)} \cdots \frac{(x-x_{j-1})}{(x_j-x_{j-1})} \frac{(x-x_{j+1})}{(x_j-x_{j+1})} \cdots \frac{(x-x_k)}{(x_j-x_k)}\$

where the product is supposed to be expanded at compile time. The \$x_m\$ are chosen to be evenly spaced on the interval \$[-1, 1]\$.

I've used boost fusion for this, mainly because it seemed (and turned out to be) helpful, but I'm not entirely sure if it's the best choice (out of the multiple TMP libraries in boost and elsewhere).

Generally, I'm new to such an extensive use of types and templates, so general advice in this regard is appreciated. And now, without further ado:

#include <iostream>

#include <boost/fusion/sequence.hpp>
#include <boost/fusion/include/accumulate.hpp>

namespace fusion = boost::fusion;


template <int Order, int NodeId>
struct Location
{
    static constexpr double coord = 2.0 * NodeId / Order - 1.0;
};


template <class BaseNode, class OtherNode>
struct Term
{
    static double F(double x)
    {
        return (x - OtherNode::coord) / (BaseNode::coord - OtherNode::coord);
    }
};


template <int Order, int BaseNodeId, int NodeCounter = Order>
struct RemainingNodes
{
    using tail = typename RemainingNodes<Order, BaseNodeId, NodeCounter - 1>::list;
    using currentNode = Location<Order, NodeCounter>;
    using full =
        typename fusion::result_of::as_vector<typename fusion::result_of::push_back<tail, currentNode>::type>::type;
    using list = typename std::conditional<BaseNodeId == NodeCounter, tail, full>::type;
};


template <int Order, int BaseNodeId>
struct RemainingNodes<Order, BaseNodeId, 0>
{
    using list = typename boost::fusion::vector<Location<Order, 0>>;
};


template <int Order>
struct RemainingNodes<Order, 0, 0>
{
    using list = typename boost::fusion::vector<>;
};


template <class BaseNode>
struct GetTerm
{
    template <class OtherNode>
    Term<BaseNode, OtherNode> operator()(OtherNode);
};


template <int Order, int BaseNodeId>
struct Terms
{
    using Nodes = typename RemainingNodes<Order, BaseNodeId>::list;
    using BaseNode = Location<Order, BaseNodeId>;
    using transformed_list = typename fusion::result_of::transform<Nodes, GetTerm<BaseNode>>::type;
    using value = typename fusion::result_of::as_vector<transformed_list>::type;
};


struct evaluate
{
    double _x;

    evaluate(double x)
        : _x(x)
    {
    }

    template <class _Term>
    double operator()(double product, const _Term&) const
    {
        return product * _Term::F(_x);
    }
};


template <int Order, int BaseNodeId>
struct LagrangePolynomial
{
    using _Terms = Terms<Order, BaseNodeId>;

    static double F(double x)
    {
        return boost::fusion::accumulate(typename _Terms::value(), 1.0, evaluate(x));
    }
};


int main()
{
    using secondOrder0 = LagrangePolynomial<2, 0>;

    std::cout << "2nd order Lagrange polynomial of 1st node, evaluated at -1, 0, 0.5, 1: ";
    std::cout << secondOrder0::F(-1.0) << " " << secondOrder0::F(0.0) << " " << secondOrder0::F(0.5) << " "
              << secondOrder0::F(1.0) << std::endl;

    return 0;
}

Answer 1

I don't want the whole function to be constexpr. I'd like to do without a for loop, but still supply the x at which to evaluate at runtime.

About constexpr

constexpr doesn't mean a function will be evaluated at compile time. It means that, given constexpr arguments, a constexpr function will be evaluated at compile-time. But that given non constexpr arguments, it will be evaluated at run-time. So, a constexpr function should give you the flexibility you want and beyond.

About Lagrange interpolation

I'll generalize a bit on your program. Let's say you don't necessarily have evenly spaced xms, just a data set of {x, y} coordinates, with the help of which you want to interpolate new coordinates. It isn't necessary to tie those two parts (coordinates generation and interpolation) too tightly, so let's not.

Let's define our types:

#include <array>

struct Point { // {x,y} coordinates
    // constexpr constructor!
    constexpr Point(double x, double y) : x(x), y(y) {}; 
    double x, y;
};

template <std::size_t N> 
using Data = std::array<Point, N>; // set of coordinates

We can then prototype our interpolation function:

template <std::size_t N>
constexpr double interpolate(const Data<N>& known_points, double n);

About compile-time loops

The last brick you need are compile-time loops. They're now trivial, thanks to "fold expressions". First, let's define a function equivalent to your term::F:

template <std::size_t N> // <=> term::F
constexpr double factor(const Data<N>& known_points, double n, int i, int j) {
// n is the x for which you want to interpolate y
// i is m in "xm"
    if (i == j) return 1;
    return (n - known_points[j].x) / (known_points[i].x - known_points[j].x);
}

You can then compute the product of all xms like this:

template <std::size_t N, std::size_t... Ns>
constexpr double factors(const Data<N>& known_points, double n, int i, std::index_sequence<Ns...>) {
    return (factor(known_points, n, i, Ns) * ...); // fold with *
}

std::index_sequence is a compile-time list of integers (more or less equivalent to boost::fusion::vectors, and a lot easier to use). Note that you can generate a std::index_sequence<0,1,..., N-1> with std::make_index_sequence<N>.

You can then complete your interpolation by summing the product of those factors by the ys:

template <std::size_t N, std::size_t... Ns>
constexpr double interpolate_impl(const Data<N>& known_points, double n, std::index_sequence<Ns...>) {
    return ( (known_points[Ns].y * factors(known_points, n, Ns, std::index_sequence<Ns...>())) + ... ); // fold with +
}

To wrap it up

Here's the program in one file to experiment with. I've added a small helper function to get a clean interface. What's important to see is that, given constexpr arguments (dataset + x whose y you need to guess), your interpolation will occur at compile-time. But you can still use it with run-time arguments -and without for loop, as requested:

#include <iostream>
#include <array>

struct Point { 
    constexpr Point(double x, double y) : x(x), y(y) {};
    double x, y;
};

template <std::size_t N>
using Data = std::array<Point, N>;


template <std::size_t N> // <=> term::F
constexpr double factor(const Data<N>& known_points, double n, int i, int j) {
    if (i == j) return 1;
    return (n - known_points[j].x) / (known_points[i].x - known_points[j].x);
}

template <std::size_t N, std::size_t... Ns>
constexpr double factors(const Data<N>& known_points, double n, int i, std::index_sequence<Ns...>) {
    return (factor(known_points, n, i, Ns)*...);
}

template <std::size_t N, std::size_t... Ns>
constexpr double interpolate_impl(const Data<N>& known_points, double n, std::index_sequence<Ns...>) {
    return ( (known_points[Ns].y * factors(known_points, n, Ns, std::index_sequence<Ns...>())) + ... );
}

template <std::size_t N>
constexpr double interpolate(const Data<N>& known_points, double n) {
    return interpolate_impl(known_points, n, std::make_index_sequence<N>());
}

int main()
{
    constexpr Data<4> kp = {Point{0.,2.}, Point{1.,3.}, Point{2.,12.}, Point{5.,147.}};
    constexpr auto r = interpolate(kp, 3.);
    static_assert(r == 35.); // completely compile-time

    double d;
    std::cin >> d;
    std::cout << interpolate(kp, d); // 35
}

Answer 2

struct evaluate ⋯

Compile-time computation is much easier today using constexpr! Rather than going for a pre-2014 TMP library, start by writing normal code and see if there is any reason why it can’t simply be declared as constexpr. Since it is a pure math function, it really should be.

I really think you can write this with no TMP whatsoever. Edit: Thanks papagaga for doing so!