constexpr class with array constructor methods

Question

I am trying to write a compile time class for multivariate polynomials (i.e. like P(X,Y,Z) = X^2 + XYZ + YZ, don't worry too much about the mathematics here):

template<int DIM, int DEGREE> class Polynomial {
  public:
    constexpr Polynomial(std::array<double,nbOfCoeffs(DIM,DEGREE)> arr): coeffs(arr) {} 


    constexpr double eval(std::array<double,DIM> x);
    constexpr operator+,-,*,/ ...
  private:
    std::array<double,nbOfCoeffs(DIM,DEGREE)> coeffs; //don't worry about "nbOfCoeffs" : it is constexpr and computes at compile time the right number of coefficients. 
}

int main () {
    Polynomial<2,2> P({{1.,1.,1.,1.,1.,1.}}); // P(X,Y) = X^2+XY+Y^2+X+Y+1

    double x = P.eval(1.);
    auto P2 = P*P;
}

So far so good no big problem to implement this. However, note that my ctor can be a bit cumbersome : how do I initialize the 3rd degree trivariate polynomial P(X,Y,Z) = XYZ ? I would have something like

Polynomial<3,3> P({{0.,0.,0.,0.,0.,1.,0.,0.,0.,0.}});

With the position of the only non-zero monomial depending on where I store it. It would be nice if I could just write :

  Polynomial<3,3> P("XYZ"); 

  Polynomial<4,3> P("X1X2X3 + X4^2"); //more general 

The idea is to create some sort of small DST to handle string representation of polynomials.

However, if I do that, once the string parsed, I don't know how to assign values to elements of my storing array: all must be done with the body of the ctor remaining empty (because I want it to be constexpr). What do you think of it ? Is it possible ? Should I change my array to some sort of recurring structure (because in this case, I think it will get really, really complex)

Answer 1

An example of how to implement Luc Danton's approach:

#include <cstddef>
#include <iostream>

namespace polynomials
{
    // it's possible to store the exponent as data member instead
    template < std::size_t t_id, std::size_t t_exponent = 1 >
    struct monomial
    {
        static constexpr std::size_t id = t_id;
        static constexpr std::size_t exponent = t_exponent;

        // it's not possible to store the coefficient
        // as non-type template parameter (floating-point..)
        double coefficient;

        explicit constexpr monomial(double p_coefficient = 1.0)
            : coefficient{ p_coefficient }
        {}

        void print() const
        {
            std::cout << coefficient << "X" << t_id << "^" << t_exponent;
        }
    };

    // create the monomial objects (like std::placeholders::_1)
    constexpr monomial<0> X;
    constexpr monomial<1> Y;
    constexpr monomial<2> Z;

    constexpr monomial<4> X0;
    constexpr monomial<5> X1;
    // ... can use macros to produce a lot of them..

    // multiply an arithmetic type (double, int, ..) with a monomial
    template < typename T_Arithmetic,
               std::size_t t_id, std::size_t t_exponent0 >
    constexpr auto operator*(T_Arithmetic c, monomial < t_id, t_exponent0 > m0)
    -> monomial < t_id, t_exponent0 >
    {
        return monomial < t_id, t_exponent0 >{c * m0.coefficient};
    }
    // the other way 'round
    template < typename T_Arithmetic,
               std::size_t t_id, std::size_t t_exponent0 >
    constexpr auto operator*(monomial < t_id, t_exponent0 > m0, T_Arithmetic c)
    -> monomial < t_id, t_exponent0 >
    {
        return c * m0;
    }

    // multiply two monomials with the same id
    template < std::size_t t_id,
               std::size_t t_exponent0, std::size_t t_exponent1 >
    constexpr auto operator*(monomial < t_id, t_exponent0 > m0,
                             monomial < t_id, t_exponent1 > m1)
    -> monomial < t_id, t_exponent0 + t_exponent1 >
    {
        return monomial<t_id, t_exponent0 + t_exponent1>
                {m0.coefficient * m1.coefficient};
    }


    // storage type for multiple different monomials
    template < typename... T_Monomials >
    struct polynomial
    {
        void print() const
        {}
    };
        template < typename T_Monomial, typename... TT_Monomials >
        struct polynomial < T_Monomial, TT_Monomials... >
            : public polynomial < TT_Monomials... >
        {
            using base = polynomial < TT_Monomials... >;

            T_Monomial m;
            constexpr polynomial(T_Monomial p, TT_Monomials... pp)
                : base(pp...)
                , m{p}
            {}

            void print() const
            {
                m.print();
                std::cout << "*";
                base::print();
            }
        };

    // multiply two monomials to get a polynomial
    template < std::size_t t_id0, std::size_t t_id1,
               std::size_t t_exponent0, std::size_t t_exponent1 >
    constexpr auto operator*( monomial < t_id0, t_exponent0 > m0,
                              monomial < t_id1, t_exponent1 > m1)
    -> polynomial < monomial<t_id0, t_exponent0>,
                    monomial<t_id1, t_exponent1> >
    {
        return {m0, m1};
    }

    // still to do (and more complicated):
    // - multiply two polynomials
    // - multiply a polynomial and a monomial
    // - addition, subtraction, division (?) etc.
}

Usage example:

int main()
{
    using namespace polynomials;

    auto p0 = 1.25*X*X;
    p0.print();
    std::cout << std::endl;

    auto p1 = p0 * 5*Y;
    p1.print();
    std::cout << std::endl;
}

Answer 2

A completely different approach to achieve the same syntax. Using member arrays vastly simplifies the development of the library:

Usage example:

int main()
{
    constexpr auto p0 = 1.25*X;
    std::cout << "p0: " << p0 << std::endl;

    constexpr auto p1 = p0*p0;
    std::cout << "p1: " << p1 << std::endl;

    constexpr auto p2 = Y*Z;  // can already multiply different monomials!!
    std::cout << "p2: " << p2 << std::endl;

    constexpr auto p3 = p1*p2;
    std::cout << "p2: " << p2 << std::endl;
}

Begin with a helper type:

#include <type_traits>
#include <iostream>

// an array type similar to `std::array`
// but with `constexpr` operators
template < typename T, std::size_t t_dim >
struct c_array
{
    T arr[t_dim];
    template < typename... TT >
    constexpr c_array(TT... pp)
        : arr{pp...}
    {}
    constexpr T operator[](std::size_t i)
    {  return arr[i]; }
    constexpr std::size_t size()
    {  return t_dim;  }
};

The monomial and polynomial types:

// the monomial type, stores a coefficient and an array of exponent
// the array index identifies a variable (X -> index 0, Y -> index 1, ..)
template < std::size_t t_numberOfVariables >
struct monomial
{
    using ExponentT = int;
    using ExponentArr = c_array < ExponentT, t_numberOfVariables >;

    double coefficient;
    ExponentArr exponents;

    // used to simplify brace-init syntax
    constexpr monomial(double c, ExponentArr e)
        : coefficient{c}
        , exponents(e)
    {}
};

// the polynomial type, stores a sum of monomials as a list (c_array)
template < std::size_t t_numberOfVariables,
           std::size_t t_numOfMonomials >
struct polynomial
{
    using MonT = monomial < t_numberOfVariables >;
    using MonArr = c_array < MonT, t_numOfMonomials >;

    MonArr monomials;

    constexpr polynomial(MonArr p)
        : monomials(p)
    {}
};

    // output / print a polynomial
    template < typename T_Char, typename T_CharTraits,
               std::size_t t_nV, std::size_t t_nP >
    std::basic_ostream<T_Char, T_CharTraits>&
    operator <<( std::basic_ostream<T_Char, T_CharTraits>& o,
                 polynomial<t_nV, t_nP> const& p )
    {
        for(std::size_t iM = 0; iM < p.monomials.size(); ++iM)
        {
            auto const& m = p.monomials[iM];

            std::cout << m.coefficient;

            for(std::size_t iExp = 0; iExp < m.exponents.size(); ++iExp)
            {
                std::cout << " * X" << iExp << "^" << m.exponents[iExp];
            }

            if( iM+1 < p.monomials.size() )
            {
                std::cout << " + ";
            }
        }

        return o;
    }

Several helpers:

// helper; construct a sequence of non-type template arguments
template < std::size_t... tt_i >
struct seq
{};

template < std::size_t t_n, std::size_t... tt_i >
struct gen_seq
    : gen_seq < t_n-1, t_n-1, tt_i...>
{};

    template < std::size_t... tt_i >
    struct gen_seq < 0, tt_i... >
        : seq < tt_i... >
    {};

// helper; compile-time max
template < typename T0, typename T1 >
constexpr auto c_max(T0 const& p0, T1 const& p1)
-> decltype( p0 > p1 ? p0 : p1 )
{
    return p0 > p1 ? p0 : p1;
}

template < typename T, typename... TT >
constexpr auto c_max(T const& p, TT const&... pp)
-> decltype( p > c_max(pp...) ? p : c_max(pp...) )
{
    return p > c_max(pp...) ? p : c_max(pp...);
}

Create the namespace-scope objects:

// helper; construct a monomial as type `polynomial`
template < std::size_t t_numberOfVariables >
constexpr polynomial<t_numberOfVariables, 1>
create_polynomial(monomial<t_numberOfVariables> m)
{
    return polynomial<t_numberOfVariables, 1>{ m };
}

template < std::size_t... tt_i >
constexpr monomial<sizeof...(tt_i) + 1>
create_monomial(double coefficient, int exponent, seq<tt_i...>)
{
    return monomial<sizeof...(tt_i) + 1>{ coefficient,
                                         {(int)(0*tt_i)..., exponent} };
}

template < std::size_t t_variableID >
constexpr polynomial<t_variableID, 1>
create_polynomial(double coefficient, int exponent)
{
    return create_polynomial<t_variableID>(
        create_monomial(coefficient, exponent, gen_seq<t_variableID-1>{}) );
}


// the namespace-scope objects
constexpr auto X  = create_monomial<1>(1.0, 1);
constexpr auto Y  = create_monomial<2>(1.0, 1);
constexpr auto Z  = create_monomial<3>(1.0, 1);
constexpr auto X0 = create_monomial<4>(1.0, 1);
// ...

Helpers for arithmetic operators on two polynomials:

// helper; expands a monomial (-> more space in array)
//         without changing its contents
template < std::size_t t_targetDim, std::size_t t_currDim,
           std::size_t... tt_curr, std::size_t... tt_exp >
constexpr monomial < t_targetDim >
expand( monomial<t_currDim> m, seq<tt_curr...>, seq<tt_exp...> )
{
    return {m.coefficient, {m.exponents[tt_curr]..., (int)(0*tt_exp)...}};
}

template < std::size_t t_targetDim, std::size_t t_currDim >
constexpr monomial < t_targetDim >
expand( monomial<t_currDim> m )
{
    using exp = std::integral_constant < std::size_t,
        (t_targetDim > t_currDim ? t_targetDim-t_currDim : 0) >;

    return expand<t_targetDim>( m, gen_seq<t_currDim>{}, gen_seq<exp{}>{} );
}

Definition of multiplication operators:

// helper for multiplication of polynomials with same array size
template < std::size_t t_dim, std::size_t... tt_i >
constexpr polynomial<t_dim>
multiply(polynomial<t_dim> p0, polynomial<t_dim> p1, seq<tt_i...>)
{
    return { p0.m[tt_i]*p1.m[tt_i]... };
}

// polynomial*polynomial, with different array size
template < std::size_t t_dim0, std::size_t t_dim1 >
constexpr polynomial < c_max(t_dim0, t_dim1) >
operator*( polynomial<t_dim0> p0, polynomial<t_dim1> p1 )
{
    using ret_dim = std::integral_constant < std::size_t,
                                             c_max(t_dim0, t_dim1) >;
    return multiply( expand<ret_dim{}>(p0), expand<ret_dim{}>(p1),
                     gen_seq<ret_dim{}>{} );
}


// helper for multiplication of monomials with same array size
template < std::size_t t_dim, std::size_t... tt_i >
constexpr monomial<t_dim>
multiply(monomial<t_dim> m0, monomial<t_dim> m1, seq<tt_i...>)
{
    return { m0.coefficient*m1.coefficient,
            {m0.exponents[tt_i]+m1.exponents[tt_i]...} };
}

// monomial*monomial, with (possibly) different array size
template < std::size_t t_dim0, std::size_t t_dim1 >
constexpr monomial < c_max(t_dim0, t_dim1) >
operator*( monomial<t_dim0> m0, monomial<t_dim1> m1 )
{
    using ret_dim = std::integral_constant < std::size_t,
                                             c_max(t_dim0, t_dim1) >;
    return multiply( expand<ret_dim{}>(m0), expand<ret_dim{}>(m1),
                     gen_seq<ret_dim{}>{} );
}


// coefficient*monomial
template < typename T_Arithmetic, std::size_t t_dim >
constexpr monomial<t_dim>
operator*(T_Arithmetic c, monomial<t_dim> m)
{
    return { c*m.coefficient, m.exponents };
}

// monomial*coefficient
template < typename T_Arithmetic, std::size_t t_dim >
constexpr monomial<t_dim>
operator*(monomial<t_dim> m, T_Arithmetic c)
{
    return { m.coefficient*c, m.exponents };
}


// helper for multiplication of coefficient*polynomial
template < typename T_Arithmetic,
           std::size_t t_nM, std::size_t t_nV,
           std::size_t... tt_i >
constexpr polynomial<t_nM, t_nV>
multiply(T_Arithmetic c, polynomial<t_nM, t_nVs> p, seq<tt_i...>)
{
    return {{c*p.monomials[tt_i]...}};
}

// helper for multiplication of polynomial*coefficient
template < typename T_Arithmetic,
           std::size_t t_nM, std::size_t t_nV,
           std::size_t... tt_i >
constexpr polynomial<t_nM, t_nV>
multiply(polynomial<t_nM, t_nV> p,
         T_Arithmetic c, seq<tt_i...>)
{
    return {{p.monomials[tt_i]*c...}};
}

// coefficient*polynomial
template < typename T_Arithmetic,
           std::size_t t_nM, std::size_t t_nV >
constexpr polynomial<t_nM, t_nV>
operator*(T_Arithmetic c, polynomial<t_nM, t_nV> p)
{
    return multiply(c, p, gen_seq<t_nM>{});
}

// polynomial*coefficient
template < typename T_Arithmetic,
           std::size_t t_nM, std::size_t t_nV >
constexpr polynomial<t_nM, t_nV>
operator*(polynomial<t_nM, t_nV> p, T_Arithmetic c)
{
    return multiply(p, c, gen_seq<t_nM>{});
}


// polynomial<N0, 1>*polynomial<N1, 1> (monomials)
template < std::size_t t_nV0,
           std::size_t t_nV1 >
constexpr polynomial< c_max(t_nV0, t_nV1), 1 >
operator*(polynomial<t_nV0, 1> p0, polynomial<t_nV1, 1> p1)
{
    return {{ p0.monomials[0]*p1.monomials[0] }};
}