Compile-time Lagrange polynomials in C++ 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 &lt;iostream&gt;

#include &lt;boost/fusion/sequence.hpp&gt;
#include &lt;boost/fusion/include/accumulate.hpp&gt;

namespace fusion = boost::fusion;


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


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


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


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


template &lt;int Order&gt;
struct RemainingNodes&lt;Order, 0, 0&gt;
{
    using list = typename boost::fusion::vector&lt;&gt;;
};


template &lt;class BaseNode&gt;
struct GetTerm
{
    template &lt;class OtherNode&gt;
    Term&lt;BaseNode, OtherNode&gt; operator()(OtherNode);
};


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


struct evaluate
{
    double _x;

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

    template &lt;class _Term&gt;
    double operator()(double product, const _Term&amp;) const
    {
        return product * _Term::F(_x);
    }
};


template &lt;int Order, int BaseNodeId&gt;
struct LagrangePolynomial
{
    using _Terms = Terms&lt;Order, BaseNodeId&gt;;

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


int main()
{
    using secondOrder0 = LagrangePolynomial&lt;2, 0&gt;;

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

    return 0;
}
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!