From b2968164e2ec74b5caa028571a10bb859f5cb5bb Mon Sep 17 00:00:00 2001 From: Theodore Turocy Date: Thu, 11 Dec 2025 12:10:11 +0000 Subject: [PATCH] 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 --- Makefile.am | 1 - src/core/array.h | 13 +- src/core/list.h | 13 +- src/solvers/enumpoly/efgpoly.cc | 4 +- src/solvers/enumpoly/poly.cc | 307 ++++++++++++++++++++++++++++++- src/solvers/enumpoly/poly.h | 234 +++++++++++++++-------- src/solvers/enumpoly/poly.imp | 291 ----------------------------- src/solvers/enumpoly/rectangle.h | 3 +- 8 files changed, 486 insertions(+), 380 deletions(-) delete mode 100644 src/solvers/enumpoly/poly.imp 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 Array { using reference = typename std::vector::reference; using const_reference = typename std::vector::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 Array { } erase(std::next(begin(), p_index - front_index())); } + template 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 T &emplace_back(Args &&...args) + { + return m_data.emplace_back(std::forward(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 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 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 BuildSequenceVariable(ProblemData &p_data, const GameSequence const std::map &var) { if (!p_sequence->action) { - return {p_data.space, 1}; + return Polynomial(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(p_data.space, var.at(p_sequence), 1); } Polynomial 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 +std::vector> Polynomial::Adder(const std::vector> &p_one, + const std::vector> &p_two) const +{ + if (p_one.empty()) { + return p_two; + } + if (p_two.empty()) { + return p_one; + } + std::vector> 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(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 +std::vector> Polynomial::Mult(const std::vector> &p_one, + const std::vector> &p_two) const +{ + std::vector> 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 &a, const Monomial &b) { return a.ExpV() < b.ExpV(); }); + std::vector> merged; + merged.reserve(result.size()); + + auto it = result.begin(); + auto end = result.end(); + + while (it != end) { + Monomial 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 Polynomial Polynomial::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 Polynomial Polynomial::PartialDerivative(int varnumber) const +{ + std::vector> 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(exponent), expv.DecrementExponent(varnumber)); + } + + Polynomial result(m_space); + result.m_terms = std::move(terms); + return result; +} + +template Polynomial Polynomial::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 Polynomial Polynomial::pow(int exponent) const +{ + if (exponent < 0) { + throw std::domain_error("Polynomial::pow: negative exponent not supported."); + } + + if (exponent == 0) { + return Polynomial(m_space, static_cast(1)); + } + if (exponent == 1) { + return *this; + } + + Polynomial base(*this); + Polynomial result(m_space, static_cast(1)); + while (exponent > 0) { + if (exponent & 1) { + result *= base; + } + exponent >>= 1; + if (exponent > 0) { + base *= base; + } + } + return result; +} + +template +Polynomial Polynomial::TranslateOfMono(const Monomial &m, + const Vector &new_origin) const +{ + Polynomial result(m_space, static_cast(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(m_space, new_origin[i]); + result *= shifted.pow(exponent); + } + result *= m.Coef(); + return result; +} + +template Polynomial Polynomial::TranslateOfPoly(const Vector &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 +Polynomial Polynomial::MonoInNewCoordinates(const Monomial &m, const Matrix &M) const +{ + Polynomial result(m_space, static_cast(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(0)) { + linearform.m_terms.emplace_back(coeff, ExponentVector(m_space, j, 1)); + } + } + result *= linearform.pow(exponent); + } + result *= m.Coef(); + return result; +} + +template Polynomial Polynomial::PolyInNewCoordinates(const Matrix &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 T Polynomial::MaximalValueOfNonlinearPart(const T &radius) const +{ + T result = static_cast(0); + + for (const auto &term : m_terms) { + if (const int deg = term.TotalDegree(); deg > 1) { + T rpow = static_cast(1); + for (int k = 0; k < deg; ++k) { + rpow *= radius; + } + result += term.Coef() * rpow; + } + } + return result; +} + +template Polynomial Polynomial::Normalize() const +{ + if (m_terms.empty()) { + return *this; + } + const T eps = std::numeric_limits::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(c > 0 ? 1 : -1)); + } + + auto it = std::max_element(m_terms.begin(), m_terms.end(), + [](const Monomial &a, const Monomial &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; template class PolynomialDerivatives; template class PolynomialSystemDerivatives; 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 m_space; Vector m_components; public: + using iterator = Vector::iterator; + using const_iterator = Vector::const_iterator; + explicit ExponentVector(std::shared_ptr 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 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 Monomial { -private: T coef; ExponentVector exps; + static T int_power(T base, int exp) + { + T result = static_cast(1); + while (exp > 0) { + if (exp & 1) { + result *= base; + } + base *= base; + exp >>= 1; + } + return result; + } + public: Monomial(std::shared_ptr 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(0)) { exps.ToZero(); } } - Monomial(const Monomial &) = default; + Monomial(const Monomial &) = default; + Monomial(Monomial &&) noexcept = default; + Monomial &operator=(Monomial &&) noexcept = default; ~Monomial() = default; // operators - Monomial &operator=(const Monomial &) = default; + Monomial &operator=(const Monomial &) = default; + + auto begin() noexcept { return exps.begin(); } + auto end() noexcept { return exps.end(); } - 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); } - Monomial operator*(const Monomial &y) const { return {coef * y.coef, exps + y.exps}; } - Monomial operator+(const Monomial &y) const { return {coef + y.coef, exps}; } - Monomial &operator+=(const Monomial &y) + auto begin() const noexcept { return exps.begin(); } + auto end() const noexcept { return exps.end(); } + + bool IsZero() const { return coef == static_cast(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 &operator*=(const T &v) + Monomial &operator*=(const T &v) { coef *= v; return *this; } - Monomial operator-() const { return {-coef, exps}; } + Monomial operator-() const { return {-coef, exps}; } std::shared_ptr 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 &vals) const + + T Evaluate(const Vector &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 Polynomial { -private: std::shared_ptr m_space; - List> m_terms; + std::vector> m_terms; // Arithmetic - List> Adder(const List> &, const List> &) const; - List> Mult(const List> &, const List> &) const; - Polynomial DivideByPolynomial(const Polynomial &den) const; + std::vector> Adder(const std::vector> &, + const std::vector> &) const; + std::vector> Mult(const std::vector> &, + const std::vector> &) const; + Polynomial DivideByPolynomial(const Polynomial &den) const; - Polynomial TranslateOfMono(const Monomial &, const Vector &) const; - Polynomial MonoInNewCoordinates(const Monomial &, const Matrix &) const; + Polynomial TranslateOfMono(const Monomial &, const Vector &) const; + Polynomial MonoInNewCoordinates(const Monomial &, const Matrix &) const; public: - Polynomial(std::shared_ptr p) : m_space(p) {} + explicit Polynomial(std::shared_ptr p) : m_space(p) {} // A constant polynomial - Polynomial(std::shared_ptr p, const T &constant) : m_space(p) + explicit Polynomial(std::shared_ptr p, const T &constant) : m_space(p) { if (constant != static_cast(0)) { m_terms.push_back(Monomial(p, constant)); } } - Polynomial(const Polynomial &) = default; + Polynomial(const Polynomial &) = default; + Polynomial(Polynomial &&) = default; + Polynomial &operator=(Polynomial &&) = default; // Constructs a polynomial x_{var_no}^exp - Polynomial(std::shared_ptr p, const int var_no, const int exp) : m_space(p) + explicit Polynomial(std::shared_ptr p, const int var_no, const int exp) + : m_space(p) { m_terms.push_back(Monomial(static_cast(1), ExponentVector(p, var_no, exp))); } // Constructs a polynomial with single monomial explicit Polynomial(const Monomial &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 &operator=(const Polynomial &) = default; - Polynomial operator-() const + Polynomial &operator=(const Polynomial &) = default; + Polynomial operator-() const { - Polynomial 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 operator-(const Polynomial &p) const + Polynomial operator-(const Polynomial &p) const { - Polynomial dif(*this); + Polynomial dif(*this); dif -= p; return dif; } - void operator-=(const Polynomial &p) + void operator-=(const Polynomial &p) { - Polynomial 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 operator+(const Polynomial &p) const + Polynomial operator+(const Polynomial &p) const { - Polynomial sum(*this); + Polynomial sum(*this); sum += p; return sum; } - Polynomial operator+(const T &v) const + Polynomial operator+(const T &v) const { - Polynomial result(*this); + Polynomial result(*this); result += v; return result; } - void operator+=(const Polynomial &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(m_space, val); } - Polynomial operator*(const Polynomial &p) const + Polynomial operator*(const Polynomial &p) const { - Polynomial prod(*this); + Polynomial prod(*this); prod *= p; return prod; } - Polynomial operator*(const T &v) const + Polynomial operator*(const T &v) const { - Polynomial result(*this); + Polynomial result(*this); result *= v; return result; } - void operator*=(const Polynomial &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 operator/(const T &val) const + Polynomial operator/(const T &val) const { if (val == static_cast(0)) { throw ZeroDivideException(); } return (*this) * (static_cast(1) / val); } - Polynomial operator/(const Polynomial &den) const { return DivideByPolynomial(den); } + Polynomial operator/(const Polynomial &den) const { return DivideByPolynomial(den); } - bool operator==(const Polynomial &p) const + bool operator==(const Polynomial &p) const { return m_space == p.m_space && m_terms == p.m_terms; } - bool operator!=(const Polynomial &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 GetSpace() const { return m_space; } - int GetDimension() const { return m_space->GetDimension(); } - int DegreeOfVar(int var_no) const + [[nodiscard]] std::shared_ptr 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 &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 &m) { return std::max(v, m.TotalDegree()); }); } - Polynomial LeadingCoefficient(int varnumber) const; + [[nodiscard]] const std::vector> &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(0); @@ -324,13 +397,12 @@ template class Polynomial { m_terms.begin(), m_terms.end(), static_cast(0), [&values](const T &v, const Monomial &m) { return v + m.Evaluate(values); }); } - Polynomial PartialDerivative(int varnumber) const; - const List> &MonomialList() const { return m_terms; } + Polynomial PartialDerivative(int varnumber) const; - Polynomial TranslateOfPoly(const Vector &) const; - Polynomial PolyInNewCoordinates(const Matrix &) const; + Polynomial TranslateOfPoly(const Vector &) const; + Polynomial PolyInNewCoordinates(const Matrix &) const; T MaximalValueOfNonlinearPart(const T &) const; - Polynomial 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 -#include "poly.h" - -namespace Gambit { - -//------------------------------------------------------------- -// Private Versions of Arithmetic Operators -//------------------------------------------------------------- - -template -List> Polynomial::Adder(const List> &One, - const List> &Two) const -{ - if (One.empty()) { - return Two; - } - if (Two.empty()) { - return One; - } - - List> 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 -List> Polynomial::Mult(const List> &One, - const List> &Two) const -{ - List> 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 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 Polynomial Polynomial::DivideByPolynomial(const Polynomial &den) const -{ - const Polynomial zero(m_space, static_cast(0)); - - if (den == zero) { - throw ZeroDivideException(); - } - // assumes exact divisibility! - - Polynomial 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 remainder = *this; - - while (remainder != zero) { - const Polynomial quot = remainder.LeadingCoefficient(last) / den.LeadingCoefficient(last); - const Polynomial 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 Polynomial Polynomial::PartialDerivative(int varnumber) const -{ - Polynomial newPoly(*this); - - for (size_t i = 1; i <= newPoly.m_terms.size(); i++) { - newPoly.m_terms[i] = - Monomial(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(0)) { - newPoly.m_terms.erase(std::next(newPoly.m_terms.begin(), j - 1)); - } - else { - j++; - } - } - - return newPoly; -} - -template Polynomial Polynomial::LeadingCoefficient(int varnumber) const -{ - Polynomial newPoly(*this); - const int degree = DegreeOfVar(varnumber); - newPoly.m_terms = List>(); - for (size_t j = 1; j <= m_terms.size(); j++) { - if (m_terms[j].ExpV()[varnumber] == degree) { - newPoly.m_terms.push_back( - Monomial(m_terms[j].Coef(), m_terms[j].ExpV().WithZeroExponent(varnumber))); - } - } - return newPoly; -} - -template -Polynomial Polynomial::TranslateOfMono(const Monomial &m, - const Vector &new_origin) const -{ - Polynomial answer(GetSpace(), static_cast(1)); - for (int i = 1; i <= GetDimension(); i++) { - if (m.ExpV()[i] > 0) { - Polynomial lt(GetSpace(), i, 1); - lt += Polynomial(GetSpace(), new_origin[i]); - for (int j = 1; j <= m.ExpV()[i]; j++) { - answer *= lt; - } - } - } - answer *= m.Coef(); - return answer; -} - -template Polynomial Polynomial::TranslateOfPoly(const Vector &new_origin) const -{ - Polynomial answer(m_space); - for (size_t i = 1; i <= this->MonomialList().size(); i++) { - answer += TranslateOfMono(this->MonomialList()[i], new_origin); - } - return answer; -} - -template -Polynomial Polynomial::MonoInNewCoordinates(const Monomial &m, const Matrix &M) const -{ - Polynomial answer(m_space, static_cast(1)); - - for (int i = 1; i <= GetDimension(); i++) { - if (m.ExpV()[i] > 0) { - Polynomial linearform(m_space, static_cast(0)); - for (int j = 1; j <= GetDimension(); j++) { - linearform += Polynomial(Monomial(M(i, j), ExponentVector(GetSpace(), j, 1))); - } - for (int k = 1; k <= m.ExpV()[i]; k++) { - answer *= linearform; - } - } - } - answer *= m.Coef(); - return answer; -} - -template Polynomial Polynomial::PolyInNewCoordinates(const Matrix &M) const -{ - Polynomial answer(m_space); - for (const auto &term : m_terms) { - answer += MonoInNewCoordinates(term, M); - } - return answer; -} - -template T Polynomial::MaximalValueOfNonlinearPart(const T &radius) const -{ - T maxcon = static_cast(0); - for (const auto &term : m_terms) { - if (term.TotalDegree() > 1) { - maxcon += term.Coef() * std::pow(radius, term.TotalDegree()); - } - } - return maxcon; -} - -template Polynomial Polynomial::Normalize() const -{ - auto maxcoeff = std::max_element( - m_terms.begin(), m_terms.end(), - [](const Monomial &a, const Monomial &b) { return a.Coef() < b.Coef(); }); - if (maxcoeff->Coef() < static_cast(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 Interval { /// A cartesian product of intervals template class Rectangle { -private: - List> sides; + Array> sides; Rectangle() = default;