Compile-time Lagrange polynomials

constexpr_lagrange_polynomials.cpp

#include <iostream>
#include <cassert>

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

namespace fusion = boost::fusion;


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


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


template <int Order, int NodeCounter = Order>
struct AllNodes
{
    using tail = typename AllNodes<Order, NodeCounter - 1>::list;
    using list = typename fusion::result_of::as_vector<
        typename fusion::result_of::push_back<tail, Location<Order, NodeCounter>>::type>::type;
};


template <int Order>
struct AllNodes<Order, 0>
{
    typedef typename boost::fusion::vector<Location<Order, 0>> list;
};


template <int Order, int BaseNode, int NodeCounter = Order>
struct RemainingNodes
{
    using tail = typename RemainingNodes<Order, BaseNode, 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<BaseNode == NodeCounter, tail, full>::type;
};


template <int Order, int BaseNode>
struct RemainingNodes<Order, BaseNode, 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 evaluateAtX
{
    double _x;

    evaluateAtX(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, evaluateAtX(x));
    }
};

struct printNodes
{
    template <typename T>
    void operator()(T const& node) const
    {
        std::cout << T::value << " ";
    }
};

struct sum
{

    template <typename T>
    double operator()(double x, const T& t) const
    {
        return x + T::value;
    }
};

struct evalAtZero
{
    template <class Term>
    void operator()(const Term&) const
    {
        std::cout << Term::F(0.0) << " ";
    }
};

int main()
{
    std::cout << "Coordinates of second node of 3rd order polynomial: " << Location<3, 1>::value << std::endl;

    std::cout << "All nodes of 2nd order polynomial: ";
    boost::fusion::for_each(AllNodes<2>::list(), printNodes());
    std::cout << std::endl;
    std::cout << "All nodes of 3nd order polynomial: ";
    boost::fusion::for_each(AllNodes<3>::list(), printNodes());
    std::cout << std::endl;
    std::cout << "All nodes of 4nd order polynomial: ";
    boost::fusion::for_each(AllNodes<4>::list(), printNodes());
    std::cout << std::endl;

    using BaseNode = Location<2, 0>;
    using OtherNode1 = Location<2, 1>;
    using OtherNode2 = Location<2, 2>;

    using firstTerm = Term<BaseNode, OtherNode1>;
    using secondTerm = Term<BaseNode, OtherNode2>;

    std::cout << "2nd order polynomial, base node is first node, evaluated at all nodes: ";
    std::cout << firstTerm::F(-1.0) * secondTerm::F(-1.0) << " ";
    std::cout << firstTerm::F(0.0) * secondTerm::F(0.0) << " ";
    std::cout << firstTerm::F(1.0) * secondTerm::F(1.0) << "\n";

    std::cout << "Nodes for 2nd order polynomial, without the 3rd node: ";
    fusion::for_each(RemainingNodes<2, 2>::list(), printNodes());
    std::cout << std::endl;
    std::cout << "Nodes for 4th order polynomial, without the 1st node: ";
    fusion::for_each(RemainingNodes<4, 0>::list(), printNodes());
    std::cout << std::endl;

    std::cout << "Sum over all nodes in 4th order polynomial, except for the 2nd: ";
    std::cout << fusion::accumulate(RemainingNodes<4, 1>::list(), 0.0, sum());
    std::cout << std::endl;

    std::cout << "Value of 2nd and 3rd term of 2nd order polynomial, evaluated at zero: ";
    fusion::for_each(Terms<2, 0>::value(), evalAtZero());
    std::cout << std::endl;

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