# Polynomial library using expression templates

In a recent toy project, I made heavy use of calculations involving polynomials in Z[x]. As an experiment, I decided to try out implementing expression templates - and it did seem to provide drastic speed-ups in my primary use case, with minimal changes to the consumer code. (In case it might be useful to see the context where I'm using this library, see: https://github.com/dschepler/groebner-zx .)

Here, my primary concern is with the expression template implementation part of the library. Some of my current thoughts (though if these particular questions are too open-ended to provide good answers, that's fair; on the other hand, if there are standard answers, in particular for the reference vs. move of subexpression objects question, then I'm all ears):

• I chose here to capture references to subexpressions throughout; so for example, auto myexpr = 2 * p + q; would immediately result in dangling references to the 2 * p subexpression. It could also be possible to move-capture intermediate subexpressions; though this would increase the implementation complexity (in particular, I would still want to avoid making copies of polynomial leaf expressions which are references to existing variables).
• One possibility I've been thinking of would be: providing a special class for subexpressions where it's easy to determine at compile time that the result of the subexpression is a monomial; and then that would allow for providing for optimized implementations of multiplication by a monomial. That would make consumer code a bit more readable: for example, instead of p -= 2 * times_x_to(q, 3); it would read p -= 2 * x_to(3) * q;. But again, at a cost of increased implementation complexity.

Things that I'm aware of but are not as much of a concern for the purposes of this review include:

• A possibility of templatizing the polynomial type so that it could also be used for polynomials over Q, over double, over mpfr arbitrary-precision reals, etc.
• Internal representation choices - in particular, whether it might be better to store the coefficients in order of increasing degree.

polynomial.h:

#pragma once

#include <gmpxx.h>
#include <initializer_list>
#include <iostream>
#include <string>
#include <vector>

using Z = mpz_class;
inline Z operator""_Z(const char* s)
{
return operator""_mpz(s);
}

// This is a library for manipulating polynomials with integer
// coefficients.  Note that because the library uses expression
// templates behind the scenes, the following general restrictions
// apply:
//
// * Use of "auto" keyword to initialize a variable with the result
//   of a calculation will generally not work as expected.  Instead,
//   explicitly set the result type to "polynomial".  e.g. instead of
//     auto p = 2 * q + r;
//   write
//     polynomial p = 2 * q + r;
//
// * Assignment expressions aliasing the destination variable will not
//   work as expected.  To work around this, use a polynomial constructor
//   to materialize the intermediate result before assignment.  e.g.
//     p += 2 * times_x_to(p + q, 3);
//   write
//     p += polynomial{2 * times_x_to(p + q, 3)};

// polynomial expression templates: each type should declare a member type
// is_polynomial_expr, and implement methods compatible with:
//    int degree_bound() const;
//      upper bound on degree of the result (does not need to be exact in
//      cases such as sum or difference of two polynomials)
//    Z coefficient(int d) const;
//      return the coefficient of x^d - where d can still be greater than
//      degree_bound()

template <typename CoeffCallable>
class polynomial_expr {
public:
using is_polynomial_expr = std::true_type;

polynomial_expr(int degree_bound, CoeffCallable&& coeff_callable)
: m_degree_bound(degree_bound)
, m_coeff_callable(std::forward<CoeffCallable>(coeff_callable))
{
}
int degree_bound() const { return m_degree_bound; }
Z coefficient(int d) const
{
return m_coeff_callable(d);
}

private:
int m_degree_bound;
typename std::decay_t<CoeffCallable> m_coeff_callable;
};

template <typename PolyExpr1, typename PolyExpr2,
typename = typename std::decay_t<PolyExpr1>::is_polynomial_expr,
typename = typename std::decay_t<PolyExpr2>::is_polynomial_expr>
auto operator+(PolyExpr1&& p, PolyExpr2&& q)
{
return polynomial_expr {
std::max(p.degree_bound(), q.degree_bound()),
[&p, &q](int d) -> Z { return p.coefficient(d) + q.coefficient(d); }
};
}

template <typename PolyExpr1, typename PolyExpr2,
typename = typename std::decay_t<PolyExpr1>::is_polynomial_expr,
typename = typename std::decay_t<PolyExpr2>::is_polynomial_expr>
auto operator-(PolyExpr1&& p, PolyExpr2&& q)
{
return polynomial_expr {
std::max(p.degree_bound(), q.degree_bound()),
[&p, &q](int d) -> Z { return p.coefficient(d) - q.coefficient(d); }
};
}

template <typename PolyExpr,
typename = typename std::decay_t<PolyExpr>::is_polynomial_expr>
auto operator-(PolyExpr&& p)
{
return polynomial_expr {
p.degree_bound(),
[&p](int d) -> Z { return -(p.coefficient(d)); }
};
}

template <typename PolyExpr,
typename = typename std::decay_t<PolyExpr>::is_polynomial_expr>
auto operator*(const Z& n, PolyExpr&& p)
{
return polynomial_expr {
n == 0 ? -1 : p.degree_bound(),
[&n, &p](int d) -> Z { return n * p.coefficient(d); }
};
}

template <typename PolyExpr,
typename = typename std::decay_t<PolyExpr>::is_polynomial_expr>
auto operator*(PolyExpr&& p, const Z& n)
{
return polynomial_expr {
n == 0 ? -1 : p.degree_bound(),
[&n, &p](int d) -> Z { return p.coefficient(d) * n; }
};
}

class polynomial;
polynomial operator*(const polynomial& p, const polynomial& q);

template <typename PolyExpr,
typename = typename std::decay_t<PolyExpr>::is_polynomial_expr>
auto times_x_to(PolyExpr&& p, int d)
{
return polynomial_expr {
p.degree_bound() < 0 ? -1 : p.degree_bound() + d,
[&p, d](int e) -> Z { return e >= d ? p.coefficient(e - d) : 0; }
};
}

template <typename PolyExpr1, typename PolyExpr2,
typename = typename std::decay_t<PolyExpr1>::is_polynomial_expr,
typename = typename std::decay_t<PolyExpr2>::is_polynomial_expr>
bool operator==(PolyExpr1&& p, PolyExpr2&& q)
{
auto d = std::max(p.degree_bound(), q.degree_bound());
for (int i = 0; i <= d; ++i)
if (p.coefficient(i) != q.coefficient(i))
return false;
return true;
}

template <typename PolyExpr1, typename PolyExpr2,
typename = typename std::decay_t<PolyExpr1>::is_polynomial_expr,
typename = typename std::decay_t<PolyExpr2>::is_polynomial_expr>
bool operator!=(PolyExpr1&& p, PolyExpr2&& q)
{
auto d = std::max(p.degree_bound(), q.degree_bound());
for (int i = 0; i <= d; ++i)
if (p.coefficient(i) != q.coefficient(i))
return true;
return false;
}

class polynomial {
public:
using is_polynomial_expr = std::true_type;

polynomial() = default;
polynomial(std::initializer_list<Z> coeffs);
explicit polynomial(std::vector<Z> coeffs);

polynomial(const polynomial&) = default;
polynomial(polynomial&&) = default;
template <typename PolyExpr,
typename = typename std::decay_t<PolyExpr>::is_polynomial_expr>
polynomial(PolyExpr&& p)
{
int d = p.degree_bound();
if (d >= 0) {
m_coeffs.reserve(d + 1);
for (; d >= 0; --d)
m_coeffs.push_back(p.coefficient(d));
normalize();
}
}

polynomial& operator=(const polynomial&) = default;
polynomial& operator=(polynomial&&) = default;
template <typename PolyExpr,
typename = typename std::decay_t<PolyExpr>::is_polynomial_expr>
polynomial& operator=(PolyExpr&& p)
{
auto deg = p.degree_bound();
m_coeffs.resize(deg + 1);
for (int d = 0; d <= deg; ++d)
m_coeffs[deg - d] = p.coefficient(d);
normalize();
return *this;
}

int degree() const { return m_coeffs.size() - 1; }
int degree_bound() const { return degree(); }
const Z& coefficient(int d) const
{
static Z static_zero = 0;
return d > degree() ? static_zero : m_coeffs[degree() - d];
}

// leading_coefficient has as a precondition that the polynomial must not be 0
{
return m_coeffs.front();
}

void negate();

template <typename PolyExpr,
typename = typename std::decay_t<PolyExpr>::is_polynomial_expr>
polynomial& operator+=(PolyExpr&& p)
{
if (p.degree_bound() > degree())
m_coeffs.insert(m_coeffs.begin(), p.degree_bound() - degree(), 0);
for (int d = p.degree_bound(); d >= 0; --d)
m_coeffs[m_coeffs.size() - d - 1] += p.coefficient(d);
normalize();
return *this;
}
template <typename PolyExpr,
typename = typename std::decay_t<PolyExpr>::is_polynomial_expr>
polynomial& operator-=(PolyExpr&& p)
{
if (p.degree_bound() > degree())
m_coeffs.insert(m_coeffs.begin(), p.degree_bound() - degree(), 0);
for (int d = p.degree_bound(); d >= 0; --d)
m_coeffs[m_coeffs.size() - d - 1] -= p.coefficient(d);
normalize();
return *this;
}
polynomial& operator*=(Z n);
polynomial& operator*=(const polynomial& p)
{
return *this = (*this * p);
}

std::string to_string() const;

friend std::ostream& operator<<(std::ostream& os, const polynomial& p)
{
return os << p.to_string();
}

private:
std::vector<Z> m_coeffs;

void normalize();
static std::string monomial_to_string(const Z& coeff, int d);
};


polynomial.cpp:

#include "polynomial.h"
#include <algorithm>

polynomial::polynomial(std::initializer_list<Z> coeffs)
: m_coeffs(coeffs)
{
normalize();
}

polynomial::polynomial(std::vector<Z> coeffs)
: m_coeffs(std::move(coeffs))
{
normalize();
}

void polynomial::normalize()
{
auto first_nonzero = find_if(m_coeffs.begin(), m_coeffs.end(),
[](const Z& coeff) { return coeff != 0; });
m_coeffs.erase(m_coeffs.begin(), first_nonzero);
}

void polynomial::negate()
{
for (auto& coeff : m_coeffs)
coeff = -coeff;
}

polynomial& polynomial::operator*=(Z n)
{
if (n == 0)
m_coeffs.clear();
else {
for (Z& coeff : m_coeffs)
coeff *= n;
}
return *this;
}

namespace polynomial_mult_details {

auto evenpart(const polynomial& p)
{
return polynomial_expr {
p.degree() / 2,
[&p](int d) -> const Z& { return p.coefficient(d * 2); }
};
}

auto oddpart(const polynomial& p)
{
return polynomial_expr {
(p.degree() - 1) / 2,
[&p](int d) -> const Z& { return p.coefficient(d * 2 + 1); }
};
}

template <typename PolyExpr1, typename PolyExpr2,
typename = typename PolyExpr1::is_polynomial_expr,
typename = typename PolyExpr2::is_polynomial_expr>
auto interleave(PolyExpr1&& p, PolyExpr2&& q)
{
return polynomial_expr {
std::max(2 * p.degree_bound(), 2 * q.degree_bound() + 1),
[&p, &q](int d) -> Z {
if (d % 2 == 0)
return p.coefficient(d / 2);
else
return q.coefficient(d / 2);
}
};
}

} // namespace polynomial_mult_details

polynomial operator*(const polynomial& p, const polynomial& q)
{
// Following the classic recursive algorithm with O(d^lg(3)) multiplications of Z values

if (p == polynomial {} || q == polynomial {})
return polynomial {};
if (p.degree() == 0)
return p.coefficient(0) * q;
if (q.degree() == 0)
return p * q.coefficient(0);

// In the following, we'll be using each coefficient of p and q
// multiple times, which is why we have designed the interface to let
// the caller materialize p and q for us.

auto pe = polynomial_mult_details::evenpart(p);
auto po = polynomial_mult_details::oddpart(p);
auto qe = polynomial_mult_details::evenpart(q);
auto qo = polynomial_mult_details::oddpart(q);

polynomial pe_qe = pe * qe;
polynomial po_qo = po * qo;
polynomial pepo_qeqo = (pe + po) * (qe + qo);

return polynomial_mult_details::interleave(
pe_qe + times_x_to(po_qo, 1),
pepo_qeqo - pe_qe - po_qo);
}

std::string polynomial::monomial_to_string(const Z& coeff, int d)
{
std::string result;
if (coeff == 1 && d > 0)
result = "";
else if (coeff == -1 && d > 0)
result = "-";
else {
result = coeff.get_str();
if (d > 0)
result += " ";
}
if (d == 1)
result += "x";
else if (d > 1) {
result += "x^";
result += std::to_string(d);
}
return result;
}

std::string polynomial::to_string() const
{
if (m_coeffs.empty())
return "0";
std::string result = monomial_to_string(m_coeffs.front(), degree());
for (int d = degree() - 1; d >= 0; --d) {
auto coeff = coefficient(d);
if (coeff > 0) {
result += " + ";
result += monomial_to_string(coeff, d);
} else if (coeff < 0) {
result += " - ";
result += monomial_to_string(-coeff, d);
}
}

return result;
}


And, in case it might be useful to see my current test battery (though I'm aware it's currently missing tests for operator+=,-=,*=): polynomial_test.cpp:

#include "polynomial.h"
#include <gtest/gtest.h>

TEST(Polynomial, Equality)
{
EXPECT_EQ((polynomial {}), (polynomial {}));
EXPECT_EQ((polynomial { 0, 0 }), (polynomial {}));
EXPECT_EQ((polynomial {}), (polynomial { 0, 0 }));
EXPECT_EQ((polynomial { 1, 2 }), (polynomial { 1, 2 }));
EXPECT_EQ((polynomial { 0, 1, 2 }), (polynomial { 1, 2 }));
EXPECT_EQ((polynomial { 1, 2 }), (polynomial { 0, 1, 2 }));
EXPECT_NE((polynomial { 1, 2 }), (polynomial { 1, 2, 0 }));
}

TEST(Polynomial, Degree)
{
EXPECT_EQ((polynomial {}.degree()), -1);
EXPECT_EQ((polynomial { 5 }.degree()), 0);
EXPECT_EQ((polynomial { 1, 2, 3 }.degree()), 2);
}

TEST(Polynomial, Coefficient)
{
EXPECT_EQ((polynomial {}.coefficient(0)), 0);
EXPECT_EQ((polynomial {}.coefficient(1)), 0);
EXPECT_EQ((polynomial { 5 }.coefficient(0)), 5);
EXPECT_EQ((polynomial { 5 }.coefficient(1)), 0);
EXPECT_EQ((polynomial { 5 }.coefficient(2)), 0);
EXPECT_EQ((polynomial { 1, 2, 3 }.coefficient(0)), 3);
EXPECT_EQ((polynomial { 1, 2, 3 }.coefficient(1)), 2);
EXPECT_EQ((polynomial { 1, 2, 3 }.coefficient(2)), 1);
EXPECT_EQ((polynomial { 1, 2, 3 }.coefficient(3)), 0);
EXPECT_EQ((polynomial { 1, 2, 3 }.coefficient(4)), 0);
}

TEST(Polynomial, Negate)
{
EXPECT_EQ((-polynomial {}), (polynomial {}));
EXPECT_EQ((-polynomial { 1, 2 }), (polynomial { -1, -2 }));
EXPECT_EQ((-polynomial { 1, -3, -2 }), (polynomial { -1, 3, 2 }));
}

TEST(Polynomial, TimesXTo)
{
EXPECT_EQ(times_x_to(polynomial {}, 5), (polynomial {}));
EXPECT_EQ(times_x_to(polynomial { 1, 2, 3 }, 0), (polynomial { 1, 2, 3 }));
EXPECT_EQ(times_x_to(polynomial { 1, 2, 3 }, 3),
(polynomial { 1, 2, 3, 0, 0, 0 }));
}

TEST(Polynomial, ScalarMult)
{
EXPECT_EQ((2 * polynomial {}), (polynomial {}));
EXPECT_EQ((0 * polynomial { 1, 2, 3 }), (polynomial {}));
EXPECT_EQ((3 * polynomial { 1, -2, 3 }), (polynomial { 3, -6, 9 }));
EXPECT_EQ((-3 * polynomial { 1, -2, 3 }), (polynomial { -3, 6, -9 }));
}

{
EXPECT_EQ((polynomial {} + polynomial {}), (polynomial {}));
EXPECT_EQ((polynomial { 1, 2, 3 } + polynomial {}), (polynomial { 1, 2, 3 }));
EXPECT_EQ((polynomial {} + polynomial { 1, 2, 3 }), (polynomial { 1, 2, 3 }));
EXPECT_EQ((polynomial { 1, 2, 3 } + polynomial { 1, 2, 3 }), (polynomial { 2, 4, 6 }));
EXPECT_EQ((polynomial { 1, 2, 3 } + polynomial { -1, -2, -3 }), (polynomial {}));
EXPECT_EQ((polynomial { 1, 2, 3, 4 } + polynomial { -1, -2, 5, 6 }),
(polynomial { 8, 10 }));
EXPECT_EQ((polynomial { 1, 2, 3, 4 } + polynomial { 5, 6 }),
(polynomial { 1, 2, 8, 10 }));
EXPECT_EQ((polynomial { 5, 6 } + polynomial { 1, 2, 3, 4 }),
(polynomial { 1, 2, 8, 10 }));
}

TEST(Polynomial, Subtract)
{
EXPECT_EQ((polynomial {} - polynomial {}), (polynomial {}));
EXPECT_EQ((polynomial { 1, 2, 3 } - polynomial {}), (polynomial { 1, 2, 3 }));
EXPECT_EQ((polynomial {} - polynomial { 1, 2, 3 }), (polynomial { -1, -2, -3 }));
EXPECT_EQ((polynomial { 1, 2, 3 } - polynomial { 1, 2, 3 }), (polynomial {}));
EXPECT_EQ((polynomial { 1, 2, 3, 4 } - polynomial { 1, 2, 5, 6 }),
(polynomial { -2, -2 }));
EXPECT_EQ((polynomial { 1, 2, 3, 4 } - polynomial { 5, 6 }),
(polynomial { 1, 2, -2, -2 }));
EXPECT_EQ((polynomial { 5, 6 } - polynomial { 1, 2, 3, 4 }),
(polynomial { -1, -2, 2, 2 }));
}

TEST(Polynomial, Multiply)
{
EXPECT_EQ((polynomial {} * polynomial {}), (polynomial {}));
EXPECT_EQ((polynomial {} * polynomial { 1, 2, 3 }), (polynomial {}));
EXPECT_EQ((polynomial { 1, 2, 3 } * polynomial {}), (polynomial {}));
EXPECT_EQ((polynomial { 2 } * polynomial { 1, 2, 3 }), (polynomial { 2, 4, 6 }));
EXPECT_EQ((polynomial { 1, 2, 3 } * polynomial { 2 }), (polynomial { 2, 4, 6 }));
EXPECT_EQ((polynomial { 1, 5 } * polynomial { 1, 6 }), (polynomial { 1, 11, 30 }));
EXPECT_EQ((polynomial { 1, 2 } * polynomial { 1, -2 }), (polynomial { 1, 0, -4 }));
EXPECT_EQ((polynomial { 1, 1, 1, 1, 1 } * polynomial { 1, 1, 1 }), (polynomial { 1, 2, 3, 3, 3, 2, 1 }));
}

TEST(Polynomial, ExprTemplates)
{
polynomial p { 2, 3, 4 };
polynomial q { 5, 6 };
EXPECT_EQ(p * q - p, p * (q - polynomial { 1 }));
EXPECT_EQ(p + q + p, 2 * p + q);
EXPECT_EQ(3 * p + 4 * times_x_to(q, 2), (polynomial { 20, 30, 9, 12 }));
EXPECT_EQ(5 * p - 2 * times_x_to(q, 1), (polynomial { 3, 20 }));
}

TEST(Polynomial, ToString)
{
EXPECT_EQ((polynomial {}.to_string()), "0");
EXPECT_EQ((polynomial { 13 }.to_string()), "13");
EXPECT_EQ((polynomial { -3 }.to_string()), "-3");
EXPECT_EQ((polynomial { 1, 0 }.to_string()), "x");
EXPECT_EQ((polynomial { 1, 4 }.to_string()), "x + 4");
EXPECT_EQ((polynomial { 1, -4 }.to_string()), "x - 4");
EXPECT_EQ((polynomial { 2, 3 }.to_string()), "2 x + 3");
EXPECT_EQ((polynomial { -2, -3 }.to_string()), "-2 x - 3");
EXPECT_EQ((polynomial { -1, 5 }.to_string()), "-x + 5");
EXPECT_EQ((polynomial { 1, 3, 2 }.to_string()), "x^2 + 3 x + 2");
EXPECT_EQ((polynomial { 1, 0, -3, 0, 0, 1, 5 }.to_string()), "x^6 - 3 x^4 + x + 5");
EXPECT_EQ((polynomial { -1, 0, 0, 0 }.to_string()), "-x^3");
EXPECT_EQ((polynomial { 1, 0, 1 }.to_string()), "x^2 + 1");
EXPECT_EQ((polynomial { 1, 0, -1 }.to_string()), "x^2 - 1");
EXPECT_EQ((polynomial { -1 }.to_string()), "-1");
}

• How dramatic were the speedups, and was that with compiler optimizations enabled or not? Oct 19, 2020 at 11:28
• In my go-to "expensive" test case, the execution time decreased from 2 seconds to 0.5 seconds. (And yes, naturally for performance comparisons I use versions with compiler optimizations enabled.) Oct 19, 2020 at 17:23

# Consider the Principle of Least Astonishment

Unless you know the implementation details of your classes, it is hard to understand why the following works:

std::cout << polynomial{1} << "\n";


But this doesn't:

std::cout << polynomial{1} + polynomial{2} << "\n";


You already explained why, but having to remember this and to wrap the expression in yet another polynomial{} is frustrating and tedious. Some users, when faced with something like this, will add polynomial{} around everything, even if it doesn't need it, and then any performance benefit will be lost.

You should have the public interface follow the principle of least astonishment, and have it work exactly like you expect regular expressions to work, as far as possible. Performance should be of secondary concern here. However, that does not mean that you shouldn't have polynomial_expr at all. It does mean however that you should try to:

• Make polynomial itself be more efficient, or
• Make polynomial_exprs implicitly cast to polynomial where needed, or
• Provide overloads for polynomial_expr that make it behave like polynomial.

For example, I can make the second statement I wrote above compile and run by adding:

template <typename CoeffCallable>
class polynomial_expr {
...
friend std::ostream &operator<<(const std::ostream &out, const polynomial_expr &expr);
};

... // after definition of class polynomial:

template <typename CoeffCallable>
std::ostream &operator<<(std::ostream &out, const polynomial_expr<CoeffCallable> &expr)
{
return out << polynomial{expr};
}


# Why are polynomial_exprs faster than polynomials?

I think the reason why returning polynomial_exprs is faster is because you avoid constructing a new polynomial with the result, which involves constructing a std::vector with the coefficients. Since that requires memory allocation, which has side-effects and might change behaviour if, for example, the global operator new is overloaded, compilers might not be able to optimize this away, even if they can see that the resulting polynomial itself is just a temporary.

But C++ users are used to this kind of performance issue; it is why we often modify containers in-place for example. So instead of writing:

polynomial p1, p2, p3 = ...;
polynomial result = p1 + p2 * p3;


If you know you no longer need to keep the original values around, you can write:

polynomial p1, p2, p3 = ...;
p1 += p2 *= p3;


While that doesn't win any beauty contests either, it is at least less surprising, and it a way a performance-conscious user can get high-performance code with just your polynomial class.

# Making polynomial go faster

As I mentioned before, the fact that polynomial uses a std::vector means it needs to do heap memory allocations. Consider using a different container for storing the coefficients that implements a small vector optimization.

Furthermore, you already mentioned making monomials a special case. I think that's a good idea, but you can maybe generalize this optimization, and make it so you only store coefficients between the highest and lowest non-zero coefficient. For example, make it so the polynomial $$\x^{102} + x^{101} + x^{100}\$$ only stores three coefficients. The implementation can be as simple as adding a member variable to class polynomial that stores the offset to the lowest non-zero coefficient. Of course, this might mean a bit more work here and there, but some operations will greatly benefit from this. For example, multiplying or dividing by a unity monomial will be trivial.

# Making it even easier to use

You already mentioned that x_to() would be a nice short way to write a monomial. Even shorter would be X(), but since that is a very short name you probably want to put it in a namespace, to avoid polluting the global namespace, and so that users can opt-in to that:

namespace polynomial_utilities {
polynomial_expr X(int d = 1) {
return {d, [&p, d](int e) -> Z { return int(e == d); }};
}
}

...

using polynomial_utilies::X;
auto expr = 3 * X(2) - 5 * X() + 4;


Alternatively you could make d a template parameter so you have to use angle brackets.

# Add a way to evaluate a polynomial at a given point

It's very nice that you can build polynomials and print them, but usually you will want to evaluate the polynomial at a given point. So I would add an operator()() to do this:

class polynomial {
...
template<typename T>
T operator()(T x) {
T result{0};
T multiplier{1};

for (auto &coeff: m_coeffs) {
result += coeff * multiplier;
multiplier *= x;
}

return result;
}
}


So you can write:

polynomial func{3, -5, 4};
std::cout << func(1.5) << "\n"; // should print: 3.25

• I suppose instead of the p1 += p2 *= p3; example, you could instead provide rvalue overloads of the operators which reuse storage (similar to what the standard library does for std::string), and then write polynomial result = std::move(p1) + std::move(p2) * std::move(p3). Right? Oct 19, 2020 at 19:49
• Though in my particular use case, the most frequent operation (as part of an implementation of modding out one polynomial in a set by another) is p -= quotient * times_x_to(q, i); where q is a const reference to an element in the set which has to remain in it; so as far as I can see, the rvalue overloads or similar wouldn't be applicable in that use case. Oct 19, 2020 at 19:53
• (And by the way, I chose to have that modding out operation be in my consumer code instead of as an operator%= in the library because in general, for the integer coefficient case, different users might want different semantics for which remainder they want - i.e. I chose to have remainder in the range [0, |leading_divisor_coeff|) whereas others might want more like standard C/C++ remainder behavior, or [-|divisor_coeff|/2, |divisor_coeff|/2], etc.) Oct 19, 2020 at 19:56
• Hmm, I'm actually a bit surprised that std::cout << polynomial{2} + polynomial{3}; doesn't work -- I would have thought it would choose the operator<<(std::ostream&, const polynomial&) overload with a conversion from the polynomial_expr<...> argument. Oct 19, 2020 at 20:21
• Oh, rvalue overloads, I hadn't thought of that. That might work too, but that requires some effort to support. As for the modulo operator, you could leave it out, but if you support it I would make it behave as much as possible as that operator does for ints, again because of least astonishment for a C++ programmer. As for operator<<(..., polynomial) not being chosen: this is because of how argument-dependent lookup works. Oct 19, 2020 at 20:44