diff --git a/Makefile.am b/Makefile.am
index 276d533b0..a1514ea45 100644
--- a/Makefile.am
+++ b/Makefile.am
@@ -404,7 +404,6 @@ gambit_enumpoly_SOURCES = \
 	src/solvers/enumpoly/rectangle.h \
 	src/solvers/enumpoly/poly.cc \
 	src/solvers/enumpoly/poly.h \
-	src/solvers/enumpoly/poly.imp \
 	src/solvers/enumpoly/polysystem.h \
     src/solvers/enumpoly/polypartial.h \
     src/solvers/enumpoly/polypartial.imp \
diff --git a/src/core/array.h b/src/core/array.h
index 57ed51936..b6e3bf1d8 100644
--- a/src/core/array.h
+++ b/src/core/array.h
@@ -51,7 +51,8 @@ template <class T> class Array {
   using reference = typename std::vector<T>::reference;
   using const_reference = typename std::vector<T>::const_reference;
 
-  explicit Array(size_t p_length = 0) : m_offset(1), m_data(p_length) {}
+  Array() : m_offset(1) {}
+  explicit Array(size_t p_length) : m_offset(1), m_data(p_length) {}
   /// Construct an array with the given index bounds.
   /// If p_high is less than p_low, the resulting array will have
   /// front_index equal to p_low, and size 0.
@@ -116,7 +117,17 @@ template <class T> class Array {
     }
     erase(std::next(begin(), p_index - front_index()));
   }
+  template <class Pred> void remove_if(Pred p)
+  {
+    auto it = std::remove_if(m_data.begin(), m_data.end(), p);
+    m_data.erase(it, m_data.end());
+  }
+
   void push_back(const T &value) { m_data.push_back(value); }
+  template <class... Args> T &emplace_back(Args &&...args)
+  {
+    return m_data.emplace_back(std::forward<Args>(args)...);
+  }
   void pop_back() { m_data.pop_back(); }
   void reserve(size_t len) { m_data.reserve(len); }
 };
diff --git a/src/core/list.h b/src/core/list.h
index be6c12e45..52dc341b1 100644
--- a/src/core/list.h
+++ b/src/core/list.h
@@ -96,7 +96,18 @@ template <class T> class List {
   /// Inserts the value at the specified position
   void insert(iterator pos, const T &value) { m_list.insert(pos, value); }
   /// Erases the element at the specified position
-  void erase(iterator pos) { m_list.erase(pos); }
+  iterator erase(iterator pos) { return m_list.erase(pos); }
+  template <class Predicate> void remove_if(Predicate pred)
+  {
+    for (auto it = begin(); it != end();) {
+      if (pred(*it)) {
+        it = erase(it);
+      }
+      else {
+        ++it;
+      }
+    }
+  }
 
   /// Removes all elements from the list container (which are destroyed),
   /// leaving the container with a size of 0.
diff --git a/src/solvers/enumpoly/efgpoly.cc b/src/solvers/enumpoly/efgpoly.cc
index e88cca353..4b642e7aa 100644
--- a/src/solvers/enumpoly/efgpoly.cc
+++ b/src/solvers/enumpoly/efgpoly.cc
@@ -60,10 +60,10 @@ Polynomial<double> BuildSequenceVariable(ProblemData &p_data, const GameSequence
                                          const std::map<GameSequence, int> &var)
 {
   if (!p_sequence->action) {
-    return {p_data.space, 1};
+    return Polynomial<double>(p_data.space, 1);
   }
   if (p_sequence->action != p_data.m_support.GetActions(p_sequence->GetInfoset()).back()) {
-    return {p_data.space, var.at(p_sequence), 1};
+    return Polynomial<double>(p_data.space, var.at(p_sequence), 1);
   }
 
   Polynomial<double> equation(p_data.space);
diff --git a/src/solvers/enumpoly/poly.cc b/src/solvers/enumpoly/poly.cc
index 80bb0e4d7..7dd79d2c4 100644
--- a/src/solvers/enumpoly/poly.cc
+++ b/src/solvers/enumpoly/poly.cc
@@ -20,11 +20,316 @@
 // Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
 //
 
-#include "poly.imp"
+#include "poly.h"
 #include "polypartial.imp"
 
 namespace Gambit {
 
+template <class T>
+std::vector<Monomial<T>> Polynomial<T>::Adder(const std::vector<Monomial<T>> &p_one,
+                                              const std::vector<Monomial<T>> &p_two) const
+{
+  if (p_one.empty()) {
+    return p_two;
+  }
+  if (p_two.empty()) {
+    return p_one;
+  }
+  std::vector<Monomial<T>> answer;
+  answer.reserve(p_one.size() + p_two.size());
+
+  auto it1 = p_one.begin();
+  auto it2 = p_two.begin();
+  auto end1 = p_one.end();
+  auto end2 = p_two.end();
+
+  while (it1 != end1 && it2 != end2) {
+    const auto &e1 = it1->ExpV();
+    const auto &e2 = it2->ExpV();
+
+    if (e1 < e2) {
+      answer.push_back(*it1++);
+    }
+    else if (e2 < e1) {
+      answer.push_back(*it2++);
+    }
+    else {
+      T merged = it1->Coef() + it2->Coef();
+      if (merged != static_cast<T>(0)) {
+        answer.emplace_back(merged, e1); // more efficient than *it1 + *it2
+      }
+      ++it1;
+      ++it2;
+    }
+  }
+
+  while (it1 != end1) {
+    answer.push_back(*it1++);
+  }
+  while (it2 != end2) {
+    answer.push_back(*it2++);
+  }
+  return answer;
+}
+
+template <class T>
+std::vector<Monomial<T>> Polynomial<T>::Mult(const std::vector<Monomial<T>> &p_one,
+                                             const std::vector<Monomial<T>> &p_two) const
+{
+  std::vector<Monomial<T>> result;
+
+  if (p_one.empty() || p_two.empty()) {
+    return result;
+  }
+
+  result.reserve(p_one.size() * p_two.size());
+  for (const auto &m1 : p_one) {
+    for (const auto &m2 : p_two) {
+      result.emplace_back(m1 * m2);
+    }
+  }
+  std::sort(result.begin(), result.end(),
+            [](const Monomial<T> &a, const Monomial<T> &b) { return a.ExpV() < b.ExpV(); });
+  std::vector<Monomial<T>> merged;
+  merged.reserve(result.size());
+
+  auto it = result.begin();
+  auto end = result.end();
+
+  while (it != end) {
+    Monomial<T> current = *it;
+    ++it;
+    while (it != end && it->ExpV() == current.ExpV()) {
+      current += *it;
+      ++it;
+    }
+    if (!current.IsZero()) {
+      merged.emplace_back(std::move(current));
+    }
+  }
+  return merged;
+}
+
+template <class T> Polynomial<T> Polynomial<T>::DivideByPolynomial(const Polynomial &den) const
+{
+  if (den.IsZero()) {
+    throw ZeroDivideException();
+  }
+
+  // assumes exact divisibility!
+  Polynomial result(m_space);
+  if (IsZero()) {
+    return result;
+  }
+
+  if (den.Degree() == 0) {
+    return *this / den.NumLeadCoeff();
+  }
+
+  int last = GetDimension();
+  while (last > 1 && den.DegreeOfVar(last) == 0) {
+    --last;
+  }
+
+  Polynomial remainder(*this);
+  while (!remainder.IsZero()) {
+    const int deg_rem = remainder.DegreeOfVar(last);
+    const int deg_den = den.DegreeOfVar(last);
+
+    const Polynomial lead_rem = remainder.LeadingCoefficient(last);
+    const Polynomial lead_den = den.LeadingCoefficient(last);
+    const Polynomial quot = lead_rem / lead_den;
+
+    // x_last^(deg_rem - deg_den)
+    const int pow = deg_rem - deg_den;
+    const Polynomial mono_pow(m_space, last, pow);
+
+    result += quot * mono_pow;
+    remainder -= quot * mono_pow * den;
+  }
+  return result;
+}
+
+template <class T> Polynomial<T> Polynomial<T>::PartialDerivative(int varnumber) const
+{
+  std::vector<Monomial<T>> terms;
+  terms.reserve(m_terms.size());
+
+  for (const auto &term : m_terms) {
+    const auto &expv = term.ExpV();
+    const int exponent = expv[varnumber];
+    if (exponent == 0) {
+      continue;
+    }
+    terms.emplace_back(term.Coef() * static_cast<T>(exponent), expv.DecrementExponent(varnumber));
+  }
+
+  Polynomial result(m_space);
+  result.m_terms = std::move(terms);
+  return result;
+}
+
+template <class T> Polynomial<T> Polynomial<T>::LeadingCoefficient(int varnumber) const
+{
+  Polynomial result(m_space);
+  result.m_terms.reserve(m_terms.size());
+
+  const int degree = DegreeOfVar(varnumber);
+  if (degree == 0) {
+    return result;
+  }
+
+  for (const auto &term : m_terms) {
+    const auto &expv = term.ExpV();
+    if (expv[varnumber] == degree) {
+      result.m_terms.emplace_back(term.Coef(), expv.WithZeroExponent(varnumber));
+    }
+  }
+  return result;
+}
+
+template <class T> Polynomial<T> Polynomial<T>::pow(int exponent) const
+{
+  if (exponent < 0) {
+    throw std::domain_error("Polynomial::pow: negative exponent not supported.");
+  }
+
+  if (exponent == 0) {
+    return Polynomial<T>(m_space, static_cast<T>(1));
+  }
+  if (exponent == 1) {
+    return *this;
+  }
+
+  Polynomial base(*this);
+  Polynomial result(m_space, static_cast<T>(1));
+  while (exponent > 0) {
+    if (exponent & 1) {
+      result *= base;
+    }
+    exponent >>= 1;
+    if (exponent > 0) {
+      base *= base;
+    }
+  }
+  return result;
+}
+
+template <class T>
+Polynomial<T> Polynomial<T>::TranslateOfMono(const Monomial<T> &m,
+                                             const Vector<T> &new_origin) const
+{
+  Polynomial result(m_space, static_cast<T>(1));
+
+  const auto &expv = m.ExpV();
+  const int dim = GetDimension();
+
+  for (int i = 1; i <= dim; ++i) {
+    const int exponent = expv[i];
+    if (exponent == 0) {
+      continue;
+    }
+    Polynomial shifted(m_space, i, 1);
+    shifted += Polynomial<T>(m_space, new_origin[i]);
+    result *= shifted.pow(exponent);
+  }
+  result *= m.Coef();
+  return result;
+}
+
+template <class T> Polynomial<T> Polynomial<T>::TranslateOfPoly(const Vector<T> &new_origin) const
+{
+  Polynomial answer(m_space);
+  answer.m_terms.reserve(m_terms.size());
+  for (const auto &mono : m_terms) {
+    answer += TranslateOfMono(mono, new_origin);
+  }
+  return answer;
+}
+
+template <class T>
+Polynomial<T> Polynomial<T>::MonoInNewCoordinates(const Monomial<T> &m, const Matrix<T> &M) const
+{
+  Polynomial result(m_space, static_cast<T>(1));
+  const auto &expv = m.ExpV();
+  const int dim = GetDimension();
+
+  int var_index = 1;
+
+  for (auto exp_it = expv.begin(); exp_it != expv.end(); ++exp_it, ++var_index) {
+    const int exponent = *exp_it;
+    if (exponent == 0) {
+      continue;
+    }
+
+    Polynomial linearform(m_space);
+    linearform.m_terms.reserve(GetDimension());
+    for (int j = 1; j <= dim; ++j) {
+      const T &coeff = M(var_index, j);
+      if (coeff != static_cast<T>(0)) {
+        linearform.m_terms.emplace_back(coeff, ExponentVector(m_space, j, 1));
+      }
+    }
+    result *= linearform.pow(exponent);
+  }
+  result *= m.Coef();
+  return result;
+}
+
+template <class T> Polynomial<T> Polynomial<T>::PolyInNewCoordinates(const Matrix<T> &M) const
+{
+  Polynomial answer(m_space);
+  answer.m_terms.reserve(m_terms.size());
+  for (const auto &term : m_terms) {
+    answer += MonoInNewCoordinates(term, M);
+  }
+  return answer;
+}
+
+template <class T> T Polynomial<T>::MaximalValueOfNonlinearPart(const T &radius) const
+{
+  T result = static_cast<T>(0);
+
+  for (const auto &term : m_terms) {
+    if (const int deg = term.TotalDegree(); deg > 1) {
+      T rpow = static_cast<T>(1);
+      for (int k = 0; k < deg; ++k) {
+        rpow *= radius;
+      }
+      result += term.Coef() * rpow;
+    }
+  }
+  return result;
+}
+
+template <class T> Polynomial<T> Polynomial<T>::Normalize() const
+{
+  if (m_terms.empty()) {
+    return *this;
+  }
+  const T eps = std::numeric_limits<T>::epsilon();
+
+  if (m_terms.size() == 1 && m_terms.front().TotalDegree() == 0) {
+    T c = m_terms.front().Coef();
+    if (std::abs(c) <= eps) {
+      return *this;
+    }
+    return Polynomial(m_space, static_cast<T>(c > 0 ? 1 : -1));
+  }
+
+  auto it = std::max_element(m_terms.begin(), m_terms.end(),
+                             [](const Monomial<T> &a, const Monomial<T> &b) {
+                               return std::abs(a.Coef()) < std::abs(b.Coef());
+                             });
+
+  const T maxc = it->Coef();
+  const T absmax = std::abs(maxc);
+  if (absmax <= eps) {
+    return *this;
+  }
+  return *this / maxc;
+}
+
 template class Polynomial<double>;
 template class PolynomialDerivatives<double>;
 template class PolynomialSystemDerivatives<double>;
diff --git a/src/solvers/enumpoly/poly.h b/src/solvers/enumpoly/poly.h
index 4bd2550e7..93cdec3a8 100644
--- a/src/solvers/enumpoly/poly.h
+++ b/src/solvers/enumpoly/poly.h
@@ -34,6 +34,7 @@ class VariableSpace {
     int number{0};
   };
 
+  VariableSpace() = delete;
   explicit VariableSpace(size_t nvars)
   {
     for (size_t i = 1; i <= nvars; i++) {
@@ -55,11 +56,13 @@ class VariableSpace {
 //
 // Exponent vectors are ordered lexicographically.
 class ExponentVector {
-private:
   std::shared_ptr<VariableSpace> m_space;
   Vector<int> m_components;
 
 public:
+  using iterator = Vector<int>::iterator;
+  using const_iterator = Vector<int>::const_iterator;
+
   explicit ExponentVector(std::shared_ptr<VariableSpace> p)
     : m_space(p), m_components(p->GetDimension())
   {
@@ -76,8 +79,14 @@ class ExponentVector {
   ~ExponentVector() = default;
   ExponentVector &operator=(const ExponentVector &) = default;
 
+  auto begin() const noexcept { return m_components.begin(); }
+  auto end() const noexcept { return m_components.end(); }
+
+  auto begin() noexcept { return m_components.begin(); }
+  auto end() noexcept { return m_components.end(); }
+
   std::shared_ptr<VariableSpace> GetSpace() const { return m_space; }
-  int operator[](int index) const { return m_components[index]; }
+  int operator[](const int index) const { return m_components[index]; }
 
   bool operator==(const ExponentVector &y) const
   {
@@ -87,6 +96,18 @@ class ExponentVector {
   {
     return m_space != y.m_space || m_components != y.m_components;
   }
+
+  bool operator<(const ExponentVector &y) const
+  {
+    if (m_space != y.m_space) {
+      throw std::logic_error("Cannot order exponent vectors from different spaces");
+    }
+    return std::lexicographical_compare(begin(), end(), y.begin(), y.end());
+  }
+  bool operator<=(const ExponentVector &y) const { return !(y < *this); }
+  bool operator>(const ExponentVector &y) const { return y < *this; }
+  bool operator>=(const ExponentVector &y) const { return !(*this < y); }
+
   ExponentVector operator+(const ExponentVector &v) const
   {
     ExponentVector tmp(*this);
@@ -100,6 +121,12 @@ class ExponentVector {
     tmp.m_components[varnumber] = 0;
     return tmp;
   }
+  ExponentVector DecrementExponent(const int varnumber) const
+  {
+    ExponentVector tmp(*this);
+    tmp.m_components[varnumber] -= 1;
+    return tmp;
+  }
 
   int GetDimension() const { return m_space->GetDimension(); }
   bool IsMultiaffine() const
@@ -109,206 +136,252 @@ class ExponentVector {
   int TotalDegree() const { return std::accumulate(m_components.begin(), m_components.end(), 0); }
 
   void ToZero() { m_components = 0; }
-
-  bool operator<(const ExponentVector &y) const
-  {
-    for (int i = 1; i <= GetDimension(); i++) {
-      if (m_components[i] < y.m_components[i]) {
-        return true;
-      }
-      if (m_components[i] > y.m_components[i]) {
-        return false;
-      }
-    }
-    return false;
-  }
-  bool operator<=(const ExponentVector &y) const { return *this < y || *this == y; }
-  bool operator>(const ExponentVector &y) const { return !(*this <= y); }
-  bool operator>=(const ExponentVector &y) const { return !(*this < y); }
 };
 
 /// A monomial of multiple variables with non-negative exponents
 template <class T> class Monomial {
-private:
   T coef;
   ExponentVector exps;
 
+  static T int_power(T base, int exp)
+  {
+    T result = static_cast<T>(1);
+    while (exp > 0) {
+      if (exp & 1) {
+        result *= base;
+      }
+      base *= base;
+      exp >>= 1;
+    }
+    return result;
+  }
+
 public:
   Monomial(std::shared_ptr<VariableSpace> p, const T &x) : coef(x), exps(p) {}
   Monomial(const T &x, const ExponentVector &e) : coef(x), exps(e)
   {
+    if (!exps.GetSpace()) {
+      throw std::logic_error("Monomial must be defined with respect to a space");
+    }
     if (x == static_cast<T>(0)) {
       exps.ToZero();
     }
   }
-  Monomial(const Monomial<T> &) = default;
+  Monomial(const Monomial &) = default;
+  Monomial(Monomial &&) noexcept = default;
+  Monomial &operator=(Monomial &&) noexcept = default;
   ~Monomial() = default;
 
   // operators
-  Monomial<T> &operator=(const Monomial<T> &) = default;
+  Monomial &operator=(const Monomial &) = default;
+
+  auto begin() noexcept { return exps.begin(); }
+  auto end() noexcept { return exps.end(); }
 
-  bool operator==(const Monomial<T> &y) const { return (coef == y.coef && exps == y.exps); }
-  bool operator!=(const Monomial<T> &y) const { return (coef != y.coef || exps != y.exps); }
-  Monomial<T> operator*(const Monomial<T> &y) const { return {coef * y.coef, exps + y.exps}; }
-  Monomial<T> operator+(const Monomial<T> &y) const { return {coef + y.coef, exps}; }
-  Monomial<T> &operator+=(const Monomial<T> &y)
+  auto begin() const noexcept { return exps.begin(); }
+  auto end() const noexcept { return exps.end(); }
+
+  bool IsZero() const { return coef == static_cast<T>(0); }
+  bool operator==(const Monomial &y) const { return (coef == y.coef && exps == y.exps); }
+  bool operator!=(const Monomial &y) const { return (coef != y.coef || exps != y.exps); }
+  bool operator<(const Monomial &y) const { return exps < y.exps; }
+  Monomial operator*(const Monomial &y) const
+  {
+    if (GetSpace() != y.GetSpace()) {
+      throw std::logic_error("Cannot multiply monomials of different spaces");
+    }
+    return {coef * y.coef, exps + y.exps};
+  }
+  Monomial operator+(const Monomial &y) const
+  {
+    if (exps != y.exps) {
+      throw std::logic_error("Attempted to add monomials with different exponents.");
+    }
+    return {coef + y.coef, exps};
+  }
+  Monomial &operator+=(const Monomial &y)
   {
+    if (exps != y.exps) {
+      throw std::logic_error("Attempted to add monomials with different exponents.");
+    }
     coef += y.coef;
     return *this;
   }
-  Monomial<T> &operator*=(const T &v)
+  Monomial &operator*=(const T &v)
   {
     coef *= v;
     return *this;
   }
-  Monomial<T> operator-() const { return {-coef, exps}; }
+  Monomial operator-() const { return {-coef, exps}; }
 
   std::shared_ptr<VariableSpace> GetSpace() const { return exps.GetSpace(); }
   const T &Coef() const { return coef; }
   int TotalDegree() const { return exps.TotalDegree(); }
   bool IsMultiaffine() const { return exps.IsMultiaffine(); }
   const ExponentVector &ExpV() const { return exps; }
-  T Evaluate(const Vector<T> &vals) const
+
+  T Evaluate(const Vector<T> &p_values) const
   {
     T answer = coef;
-    for (int i = 1; i <= exps.GetDimension(); i++) {
-      for (int j = 1; j <= exps[i]; j++) {
-        answer *= vals[i];
+
+    auto exp_it = exps.begin();
+    auto exp_end = exps.end();
+    auto val_it = p_values.begin();
+    for (; exp_it != exp_end; ++exp_it, ++val_it) {
+      if (const int e = *exp_it; e > 0) {
+        answer *= int_power(*val_it, e);
       }
     }
+
     return answer;
   }
 };
 
 // A multivariate polynomial
 template <class T> class Polynomial {
-private:
   std::shared_ptr<VariableSpace> m_space;
-  List<Monomial<T>> m_terms;
+  std::vector<Monomial<T>> m_terms;
 
   // Arithmetic
-  List<Monomial<T>> Adder(const List<Monomial<T>> &, const List<Monomial<T>> &) const;
-  List<Monomial<T>> Mult(const List<Monomial<T>> &, const List<Monomial<T>> &) const;
-  Polynomial<T> DivideByPolynomial(const Polynomial<T> &den) const;
+  std::vector<Monomial<T>> Adder(const std::vector<Monomial<T>> &,
+                                 const std::vector<Monomial<T>> &) const;
+  std::vector<Monomial<T>> Mult(const std::vector<Monomial<T>> &,
+                                const std::vector<Monomial<T>> &) const;
+  Polynomial DivideByPolynomial(const Polynomial &den) const;
 
-  Polynomial<T> TranslateOfMono(const Monomial<T> &, const Vector<T> &) const;
-  Polynomial<T> MonoInNewCoordinates(const Monomial<T> &, const Matrix<T> &) const;
+  Polynomial TranslateOfMono(const Monomial<T> &, const Vector<T> &) const;
+  Polynomial MonoInNewCoordinates(const Monomial<T> &, const Matrix<T> &) const;
 
 public:
-  Polynomial(std::shared_ptr<VariableSpace> p) : m_space(p) {}
+  explicit Polynomial(std::shared_ptr<VariableSpace> p) : m_space(p) {}
   // A constant polynomial
-  Polynomial(std::shared_ptr<VariableSpace> p, const T &constant) : m_space(p)
+  explicit Polynomial(std::shared_ptr<VariableSpace> p, const T &constant) : m_space(p)
   {
     if (constant != static_cast<T>(0)) {
       m_terms.push_back(Monomial<T>(p, constant));
     }
   }
-  Polynomial(const Polynomial<T> &) = default;
+  Polynomial(const Polynomial &) = default;
+  Polynomial(Polynomial &&) = default;
+  Polynomial &operator=(Polynomial &&) = default;
   // Constructs a polynomial x_{var_no}^exp
-  Polynomial(std::shared_ptr<VariableSpace> p, const int var_no, const int exp) : m_space(p)
+  explicit Polynomial(std::shared_ptr<VariableSpace> p, const int var_no, const int exp)
+    : m_space(p)
   {
     m_terms.push_back(Monomial<T>(static_cast<T>(1), ExponentVector(p, var_no, exp)));
   }
   // Constructs a polynomial with single monomial
   explicit Polynomial(const Monomial<T> &mono) : m_space(mono.GetSpace())
   {
+    if (!m_space) {
+      throw std::logic_error("Polynomials must be defined with respect to a space");
+    }
     m_terms.push_back(mono);
   }
   ~Polynomial() = default;
 
-  Polynomial<T> &operator=(const Polynomial<T> &) = default;
-  Polynomial<T> operator-() const
+  Polynomial &operator=(const Polynomial &) = default;
+  Polynomial operator-() const
   {
-    Polynomial<T> neg(*this);
-    for (size_t j = 1; j <= m_terms.size(); j++) {
-      neg.m_terms[j] = -m_terms[j];
+    Polynomial neg(*this);
+    for (auto &term : neg.m_terms) {
+      term = -term;
     }
     return neg;
   }
-  Polynomial<T> operator-(const Polynomial<T> &p) const
+  Polynomial operator-(const Polynomial &p) const
   {
-    Polynomial<T> dif(*this);
+    Polynomial dif(*this);
     dif -= p;
     return dif;
   }
-  void operator-=(const Polynomial<T> &p)
+  void operator-=(const Polynomial &p)
   {
-    Polynomial<T> neg = p;
-    for (size_t i = 1; i <= neg.m_terms.size(); i++) {
-      neg.m_terms[i] = -neg.m_terms[i];
+    auto neg_terms = p.m_terms;
+    for (auto &term : neg_terms) {
+      term = -term;
     }
-    m_terms = Adder(m_terms, neg.m_terms);
+    m_terms = Adder(m_terms, neg_terms);
   }
-  Polynomial<T> operator+(const Polynomial<T> &p) const
+  Polynomial operator+(const Polynomial<T> &p) const
   {
-    Polynomial<T> sum(*this);
+    Polynomial sum(*this);
     sum += p;
     return sum;
   }
-  Polynomial<T> operator+(const T &v) const
+  Polynomial operator+(const T &v) const
   {
-    Polynomial<T> result(*this);
+    Polynomial result(*this);
     result += v;
     return result;
   }
-  void operator+=(const Polynomial<T> &p) { m_terms = Adder(m_terms, p.m_terms); }
+  void operator+=(const Polynomial &p) { m_terms = Adder(m_terms, p.m_terms); }
   void operator+=(const T &val) { *this += Polynomial<T>(m_space, val); }
-  Polynomial<T> operator*(const Polynomial<T> &p) const
+  Polynomial operator*(const Polynomial &p) const
   {
-    Polynomial<T> prod(*this);
+    Polynomial prod(*this);
     prod *= p;
     return prod;
   }
-  Polynomial<T> operator*(const T &v) const
+  Polynomial operator*(const T &v) const
   {
-    Polynomial<T> result(*this);
+    Polynomial result(*this);
     result *= v;
     return result;
   }
-  void operator*=(const Polynomial<T> &p) { m_terms = Mult(m_terms, p.m_terms); }
+  void operator*=(const Polynomial &p) { m_terms = Mult(m_terms, p.m_terms); }
   void operator*=(const T &val)
   {
-    for (size_t j = 1; j <= m_terms.size(); j++) {
-      m_terms[j] *= val;
+    for (auto &term : m_terms) {
+      term *= val;
     }
   }
-  Polynomial<T> operator/(const T &val) const
+  Polynomial operator/(const T &val) const
   {
     if (val == static_cast<T>(0)) {
       throw ZeroDivideException();
     }
     return (*this) * (static_cast<T>(1) / val);
   }
-  Polynomial<T> operator/(const Polynomial<T> &den) const { return DivideByPolynomial(den); }
+  Polynomial operator/(const Polynomial &den) const { return DivideByPolynomial(den); }
 
-  bool operator==(const Polynomial<T> &p) const
+  bool operator==(const Polynomial &p) const
   {
     return m_space == p.m_space && m_terms == p.m_terms;
   }
-  bool operator!=(const Polynomial<T> &p) const
+  bool operator!=(const Polynomial &p) const
   {
     return m_space != p.m_space || m_terms != p.m_terms;
   }
 
+  Polynomial pow(int exponent) const;
+
   //-------------
   // Information
   //-------------
 
-  std::shared_ptr<VariableSpace> GetSpace() const { return m_space; }
-  int GetDimension() const { return m_space->GetDimension(); }
-  int DegreeOfVar(int var_no) const
+  [[nodiscard]] std::shared_ptr<VariableSpace> GetSpace() const noexcept { return m_space; }
+  [[nodiscard]] int GetDimension() const noexcept { return m_space->GetDimension(); }
+  [[nodiscard]] int DegreeOfVar(int var_no) const
   {
     return std::accumulate(
         m_terms.begin(), m_terms.end(), 0,
         [&var_no](int v, const Monomial<T> &m) { return std::max(v, m.ExpV()[var_no]); });
   }
-  int Degree() const
+  [[nodiscard]] int Degree() const
   {
     return std::accumulate(m_terms.begin(), m_terms.end(), 0, [](int v, const Monomial<T> &m) {
       return std::max(v, m.TotalDegree());
     });
   }
-  Polynomial<T> LeadingCoefficient(int varnumber) const;
+  [[nodiscard]] const std::vector<Monomial<T>> &GetTerms() const noexcept { return m_terms; }
+  bool IsZero() const noexcept { return m_terms.empty(); }
+  bool IsConstant() const noexcept
+  {
+    return m_terms.size() == 1 && m_terms.front().TotalDegree() == 0;
+  }
+
+  Polynomial LeadingCoefficient(int varnumber) const;
   T NumLeadCoeff() const
   {
     return (m_terms.size() == 1) ? m_terms.front().Coef() : static_cast<T>(0);
@@ -324,13 +397,12 @@ template <class T> class Polynomial {
         m_terms.begin(), m_terms.end(), static_cast<T>(0),
         [&values](const T &v, const Monomial<T> &m) { return v + m.Evaluate(values); });
   }
-  Polynomial<T> PartialDerivative(int varnumber) const;
-  const List<Monomial<T>> &MonomialList() const { return m_terms; }
+  Polynomial PartialDerivative(int varnumber) const;
 
-  Polynomial<T> TranslateOfPoly(const Vector<T> &) const;
-  Polynomial<T> PolyInNewCoordinates(const Matrix<T> &) const;
+  Polynomial TranslateOfPoly(const Vector<T> &) const;
+  Polynomial PolyInNewCoordinates(const Matrix<T> &) const;
   T MaximalValueOfNonlinearPart(const T &) const;
-  Polynomial<T> Normalize() const;
+  Polynomial Normalize() const;
 };
 
 } // end namespace Gambit
diff --git a/src/solvers/enumpoly/poly.imp b/src/solvers/enumpoly/poly.imp
deleted file mode 100644
index 9be7d8feb..000000000
--- a/src/solvers/enumpoly/poly.imp
+++ /dev/null
@@ -1,291 +0,0 @@
-//
-// This file is part of Gambit
-// Copyright (c) 1994-2025, The Gambit Project (https://www.gambit-project.org)
-//
-// FILE: src/solvers/enumpoly/poly.imp
-// Implementation of multivariate polynomial type
-//
-// This program is free software; you can redistribute it and/or modify
-// it under the terms of the GNU General Public License as published by
-// the Free Software Foundation; either version 2 of the License, or
-// (at your option) any later version.
-//
-// This program is distributed in the hope that it will be useful,
-// but WITHOUT ANY WARRANTY; without even the implied warranty of
-// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-// GNU General Public License for more details.
-//
-// You should have received a copy of the GNU General Public License
-// along with this program; if not, write to the Free Software
-// Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
-//
-
-#include <algorithm>
-#include "poly.h"
-
-namespace Gambit {
-
-//-------------------------------------------------------------
-//           Private Versions of Arithmetic Operators
-//-------------------------------------------------------------
-
-template <class T>
-List<Monomial<T>> Polynomial<T>::Adder(const List<Monomial<T>> &One,
-                                       const List<Monomial<T>> &Two) const
-{
-  if (One.empty()) {
-    return Two;
-  }
-  if (Two.empty()) {
-    return One;
-  }
-
-  List<Monomial<T>> answer;
-
-  size_t i = 1;
-  size_t j = 1;
-  while (i <= One.size() || j <= Two.size()) {
-    if (i > One.size()) {
-      answer.push_back(Two[j]);
-      j++;
-    }
-    else if (j > Two.size()) {
-      answer.push_back(One[i]);
-      i++;
-    }
-    else {
-      if (One[i].ExpV() < Two[j].ExpV()) {
-        answer.push_back(One[i]);
-        i++;
-      }
-      else if (One[i].ExpV() > Two[j].ExpV()) {
-        answer.push_back(Two[j]);
-        j++;
-      }
-      else {
-        if (One[i].Coef() + Two[j].Coef() != (T)0) {
-          answer.push_back(One[i] + Two[j]);
-        }
-        i++;
-        j++;
-      }
-    }
-  }
-
-  return answer;
-}
-
-template <class T>
-List<Monomial<T>> Polynomial<T>::Mult(const List<Monomial<T>> &One,
-                                      const List<Monomial<T>> &Two) const
-{
-  List<Monomial<T>> answer;
-
-  if (One.empty() || Two.empty()) {
-    return answer;
-  }
-
-  for (size_t i = 1; i <= One.size(); i++) {
-    for (size_t j = 1; j <= Two.size(); j++) {
-      const Monomial<T> next = One[i] * Two[j];
-
-      if (answer.empty()) {
-        answer.push_back(next);
-      }
-
-      else {
-        int bot = 1;
-        int top = answer.size();
-        if (answer[top].ExpV() < next.ExpV()) {
-          answer.push_back(next);
-        }
-        else if (answer[bot].ExpV() > next.ExpV()) {
-          answer.push_front(next);
-        }
-        else {
-          if (answer[bot].ExpV() == next.ExpV()) {
-            top = bot;
-          }
-          else if (answer[top].ExpV() == next.ExpV()) {
-            bot = top;
-          }
-
-          while (bot < top - 1) {
-            const int test = bot + (top - bot) / 2;
-            if (answer[test].ExpV() == next.ExpV()) {
-              bot = top = test;
-            }
-            else if (answer[test].ExpV() < next.ExpV()) {
-              bot = test;
-            }
-            else { // (answer[test].ExpV() > next.ExpV())
-              top = test;
-            }
-          }
-
-          if (bot == top) {
-            answer[bot] += next;
-          }
-          else {
-            answer.insert(std::next(answer.begin(), top - 1), next);
-          }
-        }
-      }
-    }
-  }
-  return answer;
-}
-
-template <class T> Polynomial<T> Polynomial<T>::DivideByPolynomial(const Polynomial<T> &den) const
-{
-  const Polynomial<T> zero(m_space, static_cast<T>(0));
-
-  if (den == zero) {
-    throw ZeroDivideException();
-  }
-  // assumes exact divisibility!
-
-  Polynomial<T> result = zero;
-
-  if (*this == zero) {
-    return result;
-  }
-  else if (den.Degree() == 0) {
-    result = *this / den.NumLeadCoeff();
-    return result;
-  }
-  else {
-    int last = GetDimension();
-    while (den.DegreeOfVar(last) == 0) {
-      last--;
-    }
-
-    Polynomial<T> remainder = *this;
-
-    while (remainder != zero) {
-      const Polynomial<T> quot = remainder.LeadingCoefficient(last) / den.LeadingCoefficient(last);
-      const Polynomial<T> power_of_last(m_space, last,
-                                        remainder.DegreeOfVar(last) - den.DegreeOfVar(last));
-      result += quot * power_of_last;
-      remainder -= quot * power_of_last * den;
-    }
-  }
-  return result;
-}
-
-template <class T> Polynomial<T> Polynomial<T>::PartialDerivative(int varnumber) const
-{
-  Polynomial<T> newPoly(*this);
-
-  for (size_t i = 1; i <= newPoly.m_terms.size(); i++) {
-    newPoly.m_terms[i] =
-        Monomial<T>(newPoly.m_terms[i].Coef() * (T)newPoly.m_terms[i].ExpV()[varnumber],
-                    newPoly.m_terms[i].ExpV().WithZeroExponent(varnumber));
-  }
-
-  size_t j = 1;
-  while (j <= newPoly.m_terms.size()) {
-    if (newPoly.m_terms[j].Coef() == static_cast<T>(0)) {
-      newPoly.m_terms.erase(std::next(newPoly.m_terms.begin(), j - 1));
-    }
-    else {
-      j++;
-    }
-  }
-
-  return newPoly;
-}
-
-template <class T> Polynomial<T> Polynomial<T>::LeadingCoefficient(int varnumber) const
-{
-  Polynomial<T> newPoly(*this);
-  const int degree = DegreeOfVar(varnumber);
-  newPoly.m_terms = List<Monomial<T>>();
-  for (size_t j = 1; j <= m_terms.size(); j++) {
-    if (m_terms[j].ExpV()[varnumber] == degree) {
-      newPoly.m_terms.push_back(
-          Monomial<T>(m_terms[j].Coef(), m_terms[j].ExpV().WithZeroExponent(varnumber)));
-    }
-  }
-  return newPoly;
-}
-
-template <class T>
-Polynomial<T> Polynomial<T>::TranslateOfMono(const Monomial<T> &m,
-                                             const Vector<T> &new_origin) const
-{
-  Polynomial<T> answer(GetSpace(), static_cast<T>(1));
-  for (int i = 1; i <= GetDimension(); i++) {
-    if (m.ExpV()[i] > 0) {
-      Polynomial<T> lt(GetSpace(), i, 1);
-      lt += Polynomial<T>(GetSpace(), new_origin[i]);
-      for (int j = 1; j <= m.ExpV()[i]; j++) {
-        answer *= lt;
-      }
-    }
-  }
-  answer *= m.Coef();
-  return answer;
-}
-
-template <class T> Polynomial<T> Polynomial<T>::TranslateOfPoly(const Vector<T> &new_origin) const
-{
-  Polynomial<T> answer(m_space);
-  for (size_t i = 1; i <= this->MonomialList().size(); i++) {
-    answer += TranslateOfMono(this->MonomialList()[i], new_origin);
-  }
-  return answer;
-}
-
-template <class T>
-Polynomial<T> Polynomial<T>::MonoInNewCoordinates(const Monomial<T> &m, const Matrix<T> &M) const
-{
-  Polynomial<T> answer(m_space, static_cast<T>(1));
-
-  for (int i = 1; i <= GetDimension(); i++) {
-    if (m.ExpV()[i] > 0) {
-      Polynomial<T> linearform(m_space, static_cast<T>(0));
-      for (int j = 1; j <= GetDimension(); j++) {
-        linearform += Polynomial<T>(Monomial<T>(M(i, j), ExponentVector(GetSpace(), j, 1)));
-      }
-      for (int k = 1; k <= m.ExpV()[i]; k++) {
-        answer *= linearform;
-      }
-    }
-  }
-  answer *= m.Coef();
-  return answer;
-}
-
-template <class T> Polynomial<T> Polynomial<T>::PolyInNewCoordinates(const Matrix<T> &M) const
-{
-  Polynomial<T> answer(m_space);
-  for (const auto &term : m_terms) {
-    answer += MonoInNewCoordinates(term, M);
-  }
-  return answer;
-}
-
-template <class T> T Polynomial<T>::MaximalValueOfNonlinearPart(const T &radius) const
-{
-  T maxcon = static_cast<T>(0);
-  for (const auto &term : m_terms) {
-    if (term.TotalDegree() > 1) {
-      maxcon += term.Coef() * std::pow(radius, term.TotalDegree());
-    }
-  }
-  return maxcon;
-}
-
-template <class T> Polynomial<T> Polynomial<T>::Normalize() const
-{
-  auto maxcoeff = std::max_element(
-      m_terms.begin(), m_terms.end(),
-      [](const Monomial<T> &a, const Monomial<T> &b) { return a.Coef() < b.Coef(); });
-  if (maxcoeff->Coef() < static_cast<T>(0.000000001)) {
-    return *this;
-  }
-  return *this / maxcoeff->Coef();
-}
-
-} // end namespace Gambit
diff --git a/src/solvers/enumpoly/rectangle.h b/src/solvers/enumpoly/rectangle.h
index 388449549..e13af9452 100644
--- a/src/solvers/enumpoly/rectangle.h
+++ b/src/solvers/enumpoly/rectangle.h
@@ -72,8 +72,7 @@ template <class T> class Interval {
 
 /// A cartesian product of intervals
 template <class T> class Rectangle {
-private:
-  List<Interval<T>> sides;
+  Array<Interval<T>> sides;
 
   Rectangle() = default;