Modernisation and efficiency improvements for polynomials.

This makes a number of updates to polynomials:
* The component monomials are now stored as a std::vector rather than a Gambit::List.
* All iteration over monomials is now iteration rather than index-based lookup (which was linear-time under List)
* Algorithm improvements have been made to addition and multiplication of
  polynomials as well as evaluation at a point.
* Revised: PartialDerivative is now valid for general polynomials (previously it was only correct for multinomials)
* Revised: Normalisation is now done if the absolute value of the largest
  coefficient exceeds a threshold (previously it was the
  value not absolute value, which was surely not intended).
* Move operations added to monomial and polynomial.

With the exception of the notes above, there should be no changes in behaviour.  Some rough benchmarking suggests this package
of changes speeds up enumpoly around 30% on some small-to-mid-size games.

LeadingCoefficient

TranslateOfMono (and adds pow)

MonoInNewCoordinates

Max value of nonlinear part

normalize

normalize - norm by absolute value

Adds IsZero, IsConstant

Use std::vector instead of Array; add explicit move to polynomial

Author: tturocy

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 \

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); }
 };

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.

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);

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>;