2

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)

Improve this question
12
  • Bear in mind that parsing strings with constexpr functions means that you will need to propagate an error somehow when constructing from an invalid run-time value. When it comes to EDSLs, have you considered an expression template approach? Commented May 19, 2013 at 19:54
  • @Luc Danton Nothing prevents you to throw a run-time exception in a constexpr. If this error in instantiated at compile time, I will get a compilation error, else it will be just as without the constexpr keyword. As for expression template, I don't know much about it, except that this is quite cumbersome to implement and that C++11 partly solves the problem more easily. I'd probably study them if I am convince I can't do it another way Commented May 19, 2013 at 20:25
  • 3
    If you feel like you can implement constexpr operators, then you can provide X, Y, Z objects and let the user go constexpr auto P = X * Y * Z;. Commented May 19, 2013 at 20:48
  • 2
    You cannot assign, but you can expand an array into a template parameter pack (-> variadic templates). Once you have the array elements as parameter pack, you can implement append, remove etc. via returning a new array that has been modified (much like in functional programming). Commented May 19, 2013 at 21:15
  • 1
    Here's an example by Xeo how to concat two arrays. Feel free to ask for different / more specific examples. Commented May 19, 2013 at 21:20

2 Answers 2

Reset to default
3

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;
}
Improve this answer
Sign up to request clarification or add additional context in comments.

1 Comment

Ok thank you very much both. I think I can definitely use this kind of API. I'll post what it looks like for complex cases (polynomial *, derivative...)
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] }};
}
Improve this answer

2 Comments

Wow ! Thank you very much. One remark : the multiply algo is wrong because you cannot multiply coeff by coeff (but I wasn't asking for the algorithm anyway). I didn't closely look yet though. Concerning the inheritance (eg gen_seq), is it purely for the empty base class optimization ?
@BérengerBerthoul Oh, you're right, there's been a design flaw as well as a mathematical flaw. I think I've fixed it now. gen_seq is just a metaprogramming tool to create an instance of seq. The latter one is another metaprogramming tool to get a sequence of numbers, which are used here as indices to expand arrays. E.g. the first multiply function. There's no need for an object of type seq, but its type provides this sequence of integers.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.