From 15b956e179ba9555e2d665561754fb1ef15bfa81 Mon Sep 17 00:00:00 2001 From: Ted Turocy Date: Thu, 12 Dec 2024 11:57:41 +0000 Subject: [PATCH] Refactor polynomial-related classes formerly in prepoly.h This reorganises and simplifies some of the preliminary classes underlying the implementation of multivariate polynomials: * Ordering of exponent vectors is now a property of vectors (or more properly the space) and not a separate concept * Functions inlined - most of the operations in these classes are simple and suitable for inlining * Rationalised and standardised namings of the concepts * Moved to gpoly.h so the polynomial implementation is now self-contained. --- Makefile.am | 8 - src/solvers/enumpoly/behavextend.cc | 87 ++- src/solvers/enumpoly/efgpoly.cc | 18 +- src/solvers/enumpoly/gpoly.cc | 20 +- src/solvers/enumpoly/gpoly.h | 289 +++++++--- src/solvers/enumpoly/gpoly.imp | 551 ++---------------- src/solvers/enumpoly/gpolylst.h | 12 +- src/solvers/enumpoly/gpolylst.imp | 17 +- src/solvers/enumpoly/ineqsolv.h | 3 +- src/solvers/enumpoly/monomial.cc | 27 - src/solvers/enumpoly/monomial.h | 66 --- src/solvers/enumpoly/monomial.imp | 134 ----- src/solvers/enumpoly/nfgpoly.cc | 26 +- src/solvers/enumpoly/poly.cc | 26 - src/solvers/enumpoly/poly.h | 150 ----- src/solvers/enumpoly/poly.imp | 862 ---------------------------- src/solvers/enumpoly/prepoly.cc | 497 ---------------- src/solvers/enumpoly/prepoly.h | 182 ------ src/solvers/enumpoly/quiksolv.h | 3 +- src/solvers/enumpoly/quiksolv.imp | 2 +- 20 files changed, 326 insertions(+), 2654 deletions(-) delete mode 100644 src/solvers/enumpoly/monomial.cc delete mode 100644 src/solvers/enumpoly/monomial.h delete mode 100644 src/solvers/enumpoly/monomial.imp delete mode 100644 src/solvers/enumpoly/poly.cc delete mode 100644 src/solvers/enumpoly/poly.h delete mode 100644 src/solvers/enumpoly/poly.imp delete mode 100644 src/solvers/enumpoly/prepoly.cc delete mode 100644 src/solvers/enumpoly/prepoly.h diff --git a/Makefile.am b/Makefile.am index 24a0abeb5..51504aac6 100644 --- a/Makefile.am +++ b/Makefile.am @@ -417,14 +417,6 @@ gambit_enumpoly_SOURCES = \ src/solvers/enumpoly/ineqsolv.h \ src/solvers/enumpoly/ineqsolv.imp \ src/solvers/enumpoly/interval.h \ - src/solvers/enumpoly/monomial.cc \ - src/solvers/enumpoly/monomial.h \ - src/solvers/enumpoly/monomial.imp \ - src/solvers/enumpoly/poly.cc \ - src/solvers/enumpoly/poly.h \ - src/solvers/enumpoly/poly.imp \ - src/solvers/enumpoly/prepoly.cc \ - src/solvers/enumpoly/prepoly.h \ src/solvers/enumpoly/quiksolv.cc \ src/solvers/enumpoly/quiksolv.h \ src/solvers/enumpoly/quiksolv.imp \ diff --git a/src/solvers/enumpoly/behavextend.cc b/src/solvers/enumpoly/behavextend.cc index 83d08d8a7..03fa78b1e 100644 --- a/src/solvers/enumpoly/behavextend.cc +++ b/src/solvers/enumpoly/behavextend.cc @@ -85,19 +85,19 @@ List DeviationInfosets(const BehaviorSupportProfile &big_supp, cons } gPolyList ActionProbsSumToOneIneqs(const MixedBehaviorProfile &p_solution, - const gSpace &BehavStratSpace, const term_order &Lex, + const VariableSpace &BehavStratSpace, const BehaviorSupportProfile &big_supp, const std::map &var_index) { - gPolyList answer(&BehavStratSpace, &Lex); + gPolyList answer(&BehavStratSpace); for (auto player : p_solution.GetGame()->GetPlayers()) { for (auto infoset : player->GetInfosets()) { if (!big_supp.HasAction(infoset)) { int index_base = var_index.at(infoset); - gPoly factor(&BehavStratSpace, 1.0, &Lex); + gPoly factor(&BehavStratSpace, 1.0); for (int k = 1; k < infoset->NumActions(); k++) { - factor -= gPoly(&BehavStratSpace, index_base + k, 1, &Lex); + factor -= gPoly(&BehavStratSpace, index_base + k, 1); } answer += factor; } @@ -177,8 +177,8 @@ std::list DeviationSupports(const BehaviorSupportProfile } bool NashNodeProbabilityPoly(const MixedBehaviorProfile &p_solution, - gPoly &node_prob, const gSpace &BehavStratSpace, - const term_order &Lex, const BehaviorSupportProfile &dsupp, + gPoly &node_prob, const VariableSpace &BehavStratSpace, + const BehaviorSupportProfile &dsupp, const std::map &var_index, GameNode tempnode, const GameInfoset &iset, const GameAction &act) { @@ -209,12 +209,12 @@ bool NashNodeProbabilityPoly(const MixedBehaviorProfile &p_solution, int initial_var_no = var_index.at(last_infoset); if (last_action->GetNumber() < last_infoset->NumActions()) { int varno = initial_var_no + last_action->GetNumber(); - node_prob *= gPoly(&BehavStratSpace, varno, 1, &Lex); + node_prob *= gPoly(&BehavStratSpace, varno, 1); } else { - gPoly factor(&BehavStratSpace, 1.0, &Lex); + gPoly factor(&BehavStratSpace, 1.0); for (int k = 1; k < last_infoset->NumActions(); k++) { - factor -= gPoly(&BehavStratSpace, initial_var_no + k, 1, &Lex); + factor -= gPoly(&BehavStratSpace, initial_var_no + k, 1); } node_prob *= factor; } @@ -225,12 +225,12 @@ bool NashNodeProbabilityPoly(const MixedBehaviorProfile &p_solution, } gPolyList NashExpectedPayoffDiffPolys(const MixedBehaviorProfile &p_solution, - const gSpace &BehavStratSpace, const term_order &Lex, + const VariableSpace &BehavStratSpace, const BehaviorSupportProfile &little_supp, const BehaviorSupportProfile &big_supp, const std::map &var_index) { - gPolyList answer(&BehavStratSpace, &Lex); + gPolyList answer(&BehavStratSpace); auto terminal_nodes = TerminalNodes(p_solution.GetGame()); @@ -249,12 +249,12 @@ gPolyList NashExpectedPayoffDiffPolys(const MixedBehaviorProfile // The utility difference between the // payoff resulting from the profile and deviation to // the strategy for pl specified by dsupp[k] - gPoly next_poly(&BehavStratSpace, &Lex); + gPoly next_poly(&BehavStratSpace); for (auto node : terminal_nodes) { - gPoly node_prob(&BehavStratSpace, 1.0, &Lex); - if (NashNodeProbabilityPoly(p_solution, node_prob, BehavStratSpace, Lex, support, - var_index, node, infoset, action)) { + gPoly node_prob(&BehavStratSpace, 1.0); + if (NashNodeProbabilityPoly(p_solution, node_prob, BehavStratSpace, support, var_index, + node, infoset, action)) { if (node->GetOutcome()) { node_prob *= static_cast(node->GetOutcome()->GetPayoff(player)); } @@ -270,15 +270,15 @@ gPolyList NashExpectedPayoffDiffPolys(const MixedBehaviorProfile } gPolyList ExtendsToNashIneqs(const MixedBehaviorProfile &p_solution, - const gSpace &BehavStratSpace, const term_order &Lex, + const VariableSpace &BehavStratSpace, const BehaviorSupportProfile &little_supp, const BehaviorSupportProfile &big_supp, const std::map &var_index) { - gPolyList answer(&BehavStratSpace, &Lex); - answer += ActionProbsSumToOneIneqs(p_solution, BehavStratSpace, Lex, big_supp, var_index); - answer += NashExpectedPayoffDiffPolys(p_solution, BehavStratSpace, Lex, little_supp, big_supp, - var_index); + gPolyList answer(&BehavStratSpace); + answer += ActionProbsSumToOneIneqs(p_solution, BehavStratSpace, big_supp, var_index); + answer += + NashExpectedPayoffDiffPolys(p_solution, BehavStratSpace, little_supp, big_supp, var_index); return answer; } @@ -306,12 +306,10 @@ bool ExtendsToNash(const MixedBehaviorProfile &p_solution, } // We establish the space - gSpace BehavStratSpace(num_vars); - ORD_PTR ptr = &lex; - term_order Lex(&BehavStratSpace, ptr); + VariableSpace BehavStratSpace(num_vars); gPolyList inequalities = - ExtendsToNashIneqs(p_solution, BehavStratSpace, Lex, little_supp, big_supp, var_index); + ExtendsToNashIneqs(p_solution, BehavStratSpace, little_supp, big_supp, var_index); // set up the rectangle of search Vector bottoms(num_vars), tops(num_vars); bottoms = 0; @@ -328,8 +326,8 @@ bool ExtendsToNash(const MixedBehaviorProfile &p_solution, namespace { bool ANFNodeProbabilityPoly(const MixedBehaviorProfile &p_solution, - gPoly &node_prob, const gSpace &BehavStratSpace, - const term_order &Lex, const BehaviorSupportProfile &big_supp, + gPoly &node_prob, const VariableSpace &BehavStratSpace, + const BehaviorSupportProfile &big_supp, const std::map &var_index, GameNode tempnode, int pl, int i, int j) { @@ -357,12 +355,12 @@ bool ANFNodeProbabilityPoly(const MixedBehaviorProfile &p_solution, int initial_var_no = var_index.at(last_infoset); if (last_action->GetNumber() < last_infoset->NumActions()) { int varno = initial_var_no + last_action->GetNumber(); - node_prob *= gPoly(&BehavStratSpace, varno, 1, &Lex); + node_prob *= gPoly(&BehavStratSpace, varno, 1); } else { - gPoly factor(&BehavStratSpace, 1.0, &Lex); + gPoly factor(&BehavStratSpace, 1.0); for (int k = 1; k < last_infoset->NumActions(); k++) { - factor -= gPoly(&BehavStratSpace, initial_var_no + k, 1, &Lex); + factor -= gPoly(&BehavStratSpace, initial_var_no + k, 1); } node_prob *= factor; } @@ -373,12 +371,12 @@ bool ANFNodeProbabilityPoly(const MixedBehaviorProfile &p_solution, } gPolyList ANFExpectedPayoffDiffPolys(const MixedBehaviorProfile &p_solution, - const gSpace &BehavStratSpace, const term_order &Lex, + const VariableSpace &BehavStratSpace, const BehaviorSupportProfile &little_supp, const BehaviorSupportProfile &big_supp, const std::map &var_index) { - gPolyList answer(&BehavStratSpace, &Lex); + gPolyList answer(&BehavStratSpace); auto terminal_nodes = TerminalNodes(p_solution.GetGame()); @@ -394,12 +392,12 @@ gPolyList ANFExpectedPayoffDiffPolys(const MixedBehaviorProfile // This will be the utility difference between the // payoff resulting from the profile and deviation to // action j - gPoly next_poly(&BehavStratSpace, &Lex); + gPoly next_poly(&BehavStratSpace); for (auto terminal : terminal_nodes) { - gPoly node_prob(&BehavStratSpace, 1.0, &Lex); - if (ANFNodeProbabilityPoly(p_solution, node_prob, BehavStratSpace, Lex, big_supp, - var_index, terminal, player->GetNumber(), - infoset->GetNumber(), action->GetNumber())) { + gPoly node_prob(&BehavStratSpace, 1.0); + if (ANFNodeProbabilityPoly(p_solution, node_prob, BehavStratSpace, big_supp, var_index, + terminal, player->GetNumber(), infoset->GetNumber(), + action->GetNumber())) { if (terminal->GetOutcome()) { node_prob *= static_cast(terminal->GetOutcome()->GetPayoff(player)); } @@ -414,15 +412,15 @@ gPolyList ANFExpectedPayoffDiffPolys(const MixedBehaviorProfile } gPolyList ExtendsToANFNashIneqs(const MixedBehaviorProfile &p_solution, - const gSpace &BehavStratSpace, const term_order &Lex, + const VariableSpace &BehavStratSpace, const BehaviorSupportProfile &little_supp, const BehaviorSupportProfile &big_supp, const std::map &var_index) { - gPolyList answer(&BehavStratSpace, &Lex); - answer += ActionProbsSumToOneIneqs(p_solution, BehavStratSpace, Lex, big_supp, var_index); - answer += ANFExpectedPayoffDiffPolys(p_solution, BehavStratSpace, Lex, little_supp, big_supp, - var_index); + gPolyList answer(&BehavStratSpace); + answer += ActionProbsSumToOneIneqs(p_solution, BehavStratSpace, big_supp, var_index); + answer += + ANFExpectedPayoffDiffPolys(p_solution, BehavStratSpace, little_supp, big_supp, var_index); return answer; } @@ -448,13 +446,10 @@ bool ExtendsToAgentNash(const MixedBehaviorProfile &p_solution, } // We establish the space - gSpace BehavStratSpace(num_vars); - ORD_PTR ptr = &lex; - term_order Lex(&BehavStratSpace, ptr); - + VariableSpace BehavStratSpace(num_vars); num_vars = BehavStratSpace.Dmnsn(); gPolyList inequalities = - ExtendsToANFNashIneqs(p_solution, BehavStratSpace, Lex, little_supp, big_supp, var_index); + ExtendsToANFNashIneqs(p_solution, BehavStratSpace, little_supp, big_supp, var_index); // set up the rectangle of search Vector bottoms(num_vars), tops(num_vars); diff --git a/src/solvers/enumpoly/efgpoly.cc b/src/solvers/enumpoly/efgpoly.cc index 03d5537e5..c4bf5200d 100644 --- a/src/solvers/enumpoly/efgpoly.cc +++ b/src/solvers/enumpoly/efgpoly.cc @@ -53,8 +53,7 @@ namespace { class ProblemData { public: GameSequenceForm sfg; - gSpace Space; - term_order Lex; + VariableSpace Space; std::map var; std::map> variables; @@ -65,13 +64,13 @@ gPoly BuildSequenceVariable(ProblemData &p_data, const GameSequence &p_s const std::map &var) { if (!p_sequence->action) { - return {&p_data.Space, 1, &p_data.Lex}; + return {&p_data.Space, 1}; } if (p_sequence->action != p_data.sfg.GetSupport().GetActions(p_sequence->GetInfoset()).back()) { - return {&p_data.Space, var.at(p_sequence), 1, &p_data.Lex}; + return {&p_data.Space, var.at(p_sequence), 1}; } - gPoly equation(&p_data.Space, &p_data.Lex); + gPoly equation(&p_data.Space); for (auto seq : p_data.sfg.GetSequences(p_sequence->player)) { if (seq == p_sequence) { continue; @@ -86,8 +85,7 @@ gPoly BuildSequenceVariable(ProblemData &p_data, const GameSequence &p_s ProblemData::ProblemData(const BehaviorSupportProfile &p_support) : sfg(p_support), - Space(sfg.GetSequences().size() - sfg.GetInfosets().size() - sfg.GetPlayers().size()), - Lex(&Space, lex) + Space(sfg.GetSequences().size() - sfg.GetInfosets().size() - sfg.GetPlayers().size()) { for (auto sequence : sfg.GetSequences()) { if (sequence->action && @@ -103,12 +101,12 @@ ProblemData::ProblemData(const BehaviorSupportProfile &p_support) gPoly GetPayoff(ProblemData &p_data, const GamePlayer &p_player) { - gPoly equation(&p_data.Space, &p_data.Lex); + gPoly equation(&p_data.Space); for (auto profile : p_data.sfg.GetContingencies()) { auto pay = p_data.sfg.GetPayoff(profile, p_player); if (pay != Rational(0)) { - gPoly term(&p_data.Space, double(pay), &p_data.Lex); + gPoly term(&p_data.Space, double(pay)); for (auto player : p_data.sfg.GetPlayers()) { term *= p_data.variables.at(profile[player]); } @@ -160,7 +158,7 @@ std::list> SolveSupport(const BehaviorSupportProfil bool &p_isSingular, int p_stopAfter) { ProblemData data(p_support); - gPolyList equations(&data.Space, &data.Lex); + gPolyList equations(&data.Space); IndifferenceEquations(data, equations); LastActionProbPositiveInequalities(data, equations); diff --git a/src/solvers/enumpoly/gpoly.cc b/src/solvers/enumpoly/gpoly.cc index 41c90646b..e7a8374e1 100644 --- a/src/solvers/enumpoly/gpoly.cc +++ b/src/solvers/enumpoly/gpoly.cc @@ -23,28 +23,10 @@ #include "gpoly.imp" #include "gambit.h" -template <> double gPoly::String_Coeff(double nega) -{ - std::string Coeff; - while ((charc >= '0' && charc <= '9') || charc == '.') { - Coeff += charc; - charnum++; - GetChar(); - } - if (Coeff.empty()) { - return nega; - } - else { - return (nega * strtod(Coeff.c_str(), nullptr)); - } -} - template class gPoly; template gPoly operator+(const gPoly &poly, const double &val); template gPoly operator*(const double &val, const gPoly &poly); template gPoly operator*(const gPoly &poly, const double &val); -template gPoly TogDouble(const gPoly &); +template gPoly ToDouble(const gPoly &); template gPoly NormalizationOfPoly(const gPoly &); - -template std::string &operator<<(std::string &, const gPoly &); diff --git a/src/solvers/enumpoly/gpoly.h b/src/solvers/enumpoly/gpoly.h index 56bb16324..28dd2fab7 100644 --- a/src/solvers/enumpoly/gpoly.h +++ b/src/solvers/enumpoly/gpoly.h @@ -25,38 +25,215 @@ #include "gambit.h" #include "core/sqmatrix.h" -#include "monomial.h" -#include "poly.h" -// These classes are used to store and mathematically manipulate polynomials. +using namespace Gambit; + +class VariableSpace { +public: + struct Variable { + std::string name; + int number; + }; + + explicit VariableSpace(size_t nvars) + { + for (size_t i = 1; i <= nvars; i++) { + m_variables.push_back({"n" + std::to_string(i), static_cast(i)}); + } + } + VariableSpace(const VariableSpace &) = delete; + ~VariableSpace() = default; + VariableSpace &operator=(const VariableSpace &) = delete; + const Variable &operator[](const int index) const { return m_variables[index]; } + int Dmnsn() const { return m_variables.size(); } + +private: + Array m_variables; +}; + +// An exponent vector is a vector of integers representing the exponents on variables in +// a space in a monomial. +// +// Exponent vectors are ordered lexicographically. +class ExponentVector { +private: + const VariableSpace *m_space; + Vector m_components; + +public: + explicit ExponentVector(const VariableSpace *p) : m_space(p), m_components(p->Dmnsn()) + { + m_components = 0; + } + // Construct x_i^j + ExponentVector(const VariableSpace *p, const int i, const int j) + : m_space(p), m_components(p->Dmnsn()) + { + m_components = 0; + m_components[i] = j; + } + ExponentVector(const VariableSpace *space, const Vector &exps) + : m_space(space), m_components(exps) + { + } + ExponentVector(const ExponentVector &) = default; + ~ExponentVector() = default; + + // Operators + ExponentVector &operator=(const ExponentVector &) = default; + + int operator[](int index) const { return m_components[index]; } + + bool operator==(const ExponentVector &y) const + { + return m_space == y.m_space && m_components == y.m_components; + } + bool operator!=(const ExponentVector &y) const + { + return m_space != y.m_space || m_components != y.m_components; + } + ExponentVector operator+(const ExponentVector &v) const + { + ExponentVector tmp(*this); + tmp.m_components += v.m_components; + return tmp; + } + + // Other operations + ExponentVector WithZeroExponent(const int varnumber) const + { + ExponentVector tmp(*this); + tmp.m_components[varnumber] = 0; + return tmp; + } + ExponentVector WithDecrementedExponent(const int varnumber) const + { + ExponentVector tmp(*this); + tmp.m_components[varnumber]--; + return tmp; + } + + // Information + int Dmnsn() const { return m_space->Dmnsn(); } + bool IsConstant() const + { + for (int i = 1; i <= Dmnsn(); i++) { + if ((*this)[i] > 0) { + return false; + } + } + return true; + } + bool IsMultiaffine() const + { + for (int i = 1; i <= Dmnsn(); i++) { + if ((*this)[i] > 1) { + return false; + } + } + return true; + } + int TotalDegree() const + { + int exp_sum = 0; + for (int i = 1; i <= Dmnsn(); i++) { + exp_sum += (*this)[i]; + } + return exp_sum; + } + + // Manipulation + void ToZero() + { + for (int i = 1; i <= Dmnsn(); i++) { + m_components[i] = 0; + } + } + + bool operator<(const ExponentVector &y) const + { + for (int i = 1; i <= Dmnsn(); 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); } +}; -// **NOTE** -// Every type T to be used needs a procedure to convert a gText coefficient -// to the type T for the gText SOP input form and a procedure to convert -// the coefficient into a gText for the SOP output form. +/// A monomial of multiple variables with non-negative exponents +template class gMono { +private: + T coef; + ExponentVector exps; -// ******************* -// gPoly declaration -// ******************* +public: + // constructors + gMono(const VariableSpace *p, const T &x) : coef(x), exps(p) {} + gMono(const T &x, const ExponentVector &e) : coef(x), exps(e) + { + if (x == static_cast(0)) { + exps.ToZero(); + } + } + gMono(const gMono &) = default; + ~gMono() = default; + + // operators + gMono &operator=(const gMono &) = default; + + bool operator==(const gMono &y) const { return (coef == y.coef && exps == y.exps); } + bool operator!=(const gMono &y) const { return (coef != y.coef || exps != y.exps); } + gMono operator*(const gMono &y) const { return {coef * y.coef, exps + y.exps}; } + gMono operator+(const gMono &y) const { return {coef + y.coef, exps}; } + gMono &operator+=(const gMono &y) + { + coef += y.coef; + return *this; + } + gMono &operator*=(const T &v) + { + coef *= v; + return *this; + } + gMono operator-() const { return {-coef, exps}; } + + // information + const T &Coef() const { return coef; } + int Dmnsn() const { return exps.Dmnsn(); } + int TotalDegree() const { return exps.TotalDegree(); } + bool IsConstant() const { return exps.IsConstant(); } + bool IsMultiaffine() const { return exps.IsMultiaffine(); } + const ExponentVector &ExpV() const { return exps; } + T Evaluate(const Vector &vals) const + { + T answer = coef; + for (int i = 1; i <= exps.Dmnsn(); i++) { + for (int j = 1; j <= exps[i]; j++) { + answer *= vals[i]; + } + } + return answer; + } +}; +// These classes are used to store and mathematically manipulate polynomials. template class gPoly { private: - const gSpace *Space; // pointer to variable Space of space - const term_order *Order; + const VariableSpace *Space; // pointer to variable Space of space Gambit::List> Terms; // alternative implementation - // used for gText parsing; - unsigned int charnum; - char charc; - std::string TheString; - //---------------------- // some private members //---------------------- - // Information - exp_vect OrderMaxMonomialDivisibleBy(const term_order &order, const exp_vect &expv); // Arithmetic Gambit::List> Adder(const Gambit::List> &, const Gambit::List> &) const; @@ -68,45 +245,23 @@ template class gPoly { gPoly TranslateOfMono(const gMono &, const Gambit::Vector &) const; gPoly MonoInNewCoordinates(const gMono &, const Gambit::SquareMatrix &) const; - //----------------------------------------------- - // Going back and forth from std::strings to gPoly's - //----------------------------------------------- - - // std::string input parser functions - void String_Term(T nega); - T String_Coeff(T nega); - int String_GetPow(); - void String_VarAndPow(Gambit::Array &PowArray); - void GetChar(); - // Is the string a valid polynomial? - bool Check_String(const std::string &Hold); - - //---------------------- - // private friends - //---------------------- - - // friend gPoly operator*<>(const gPoly &poly, const T val); - // friend gPoly operator*(const T val, const gPoly &poly); - public: //--------------------------- // Construction, destruction: //--------------------------- // Null gPoly constructor - gPoly(const gSpace *, const term_order *); - // Constructs a gPoly equal to the SOP representation in the std::string - gPoly(const gSpace *, const std::string &, const term_order *); + gPoly(const VariableSpace *); // Constructs a constant gPoly - gPoly(const gSpace *, const T &, const term_order *); + gPoly(const VariableSpace *, const T &); // Constructs a gPoly equal to another; gPoly(const gPoly &); // Constructs a gPoly that is x_{var_no}^exp; - gPoly(const gSpace *p, int var_no, int exp, const term_order *); + gPoly(const VariableSpace *p, int var_no, int exp); // Constructs a gPoly that is the monomial coeff*vars^exps; - gPoly(const gSpace *p, exp_vect exps, T coeff, const term_order *); + gPoly(const VariableSpace *p, ExponentVector exps, T coeff); // Constructs a gPoly with single monomial - gPoly(const gSpace *p, const gMono &, const term_order *); + gPoly(const VariableSpace *p, const gMono &); ~gPoly() = default; @@ -115,7 +270,6 @@ template class gPoly { //---------- gPoly &operator=(const gPoly &); - gPoly &operator=(const std::string &); // Set polynomial equal to the SOP form in the string gPoly operator-() const; gPoly operator-(const gPoly &) const; @@ -136,60 +290,31 @@ template class gPoly { // Information: //------------- - const gSpace *GetSpace() const; - const term_order *GetOrder() const; + const VariableSpace *GetSpace() const; int Dmnsn() const; - bool IsZero() const; int DegreeOfVar(int var_no) const; int Degree() const; - T GetCoef(const Gambit::Array &Powers) const; - T GetCoef(const exp_vect &Powers) const; + T GetCoef(const Gambit::Vector &Powers) const; + T GetCoef(const ExponentVector &Powers) const; gPoly LeadingCoefficient(int varnumber) const; T NumLeadCoeff() const; // deg == 0 bool IsConstant() const; bool IsMultiaffine() const; - int UniqueActiveVariable() const; - // returns 0 if constant, -1 if truly multivariate - polynomial UnivariateEquivalent(int activar) const; // assumes UniqueActiveVariable() is true - T Evaluate(const Gambit::Array &values) const; - gPoly EvaluateOneVar(int varnumber, T val) const; + T Evaluate(const Gambit::Vector &values) const; gPoly PartialDerivative(int varnumber) const; - int No_Monomials() const; - Gambit::List ExponentVectors() const; Gambit::List> MonomialList() const; gPoly TranslateOfPoly(const Gambit::Vector &) const; gPoly PolyInNewCoordinates(const Gambit::SquareMatrix &) const; T MaximalValueOfNonlinearPart(const T &) const; - - //-------------------- - // Term Order Concepts - //-------------------- - - exp_vect LeadingPowerProduct(const term_order &) const; - T LeadingCoefficient(const term_order &) const; - gPoly LeadingTerm(const term_order &) const; - void ToMonic(const term_order &); - void ReduceByDivisionAtExpV(const term_order &, const gPoly &, const exp_vect &); - void ReduceByRepeatedDivision(const term_order &, const gPoly &); - gPoly S_Polynomial(const term_order &, const gPoly &) const; - - //--------------- - // Printing Stuff - //--------------- - - // Print polynomial in SOP form - void Output(std::string &) const; }; -template std::string &operator<<(std::string &, const gPoly &); - //------------- // Conversion: //------------- -template gPoly TogDouble(const gPoly &); +template gPoly ToDouble(const gPoly &); template gPoly NormalizationOfPoly(const gPoly &); // global multiply by scalar operators @@ -200,6 +325,4 @@ template gPoly operator*(const gPoly &poly, const T &val); template gPoly operator+(const T &val, const gPoly &poly); template gPoly operator+(const gPoly &poly, const T &val); -template std::string ToText(const gPoly &p); - #endif // # GPOLY_H diff --git a/src/solvers/enumpoly/gpoly.imp b/src/solvers/enumpoly/gpoly.imp index 6d0ab3ca5..c402cfb38 100644 --- a/src/solvers/enumpoly/gpoly.imp +++ b/src/solvers/enumpoly/gpoly.imp @@ -20,7 +20,6 @@ // Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. // -#include // for abs() #include // for std::max() #include "gpoly.h" #include "gambit.h" @@ -34,13 +33,13 @@ //--------------------------- template -gPoly::gPoly(const gSpace *p, const term_order *o) : Space(p), Order(o), Terms() +gPoly::gPoly(const VariableSpace *p) : Space(p) { } template -gPoly::gPoly(const gSpace *p, const T &constant, const term_order *o) - : Space(p), Order(o), Terms() +gPoly::gPoly(const VariableSpace *p, const T &constant) + : Space(p) { if (constant != (T)0) { Terms.push_back(gMono(p, constant)); @@ -48,35 +47,28 @@ gPoly::gPoly(const gSpace *p, const T &constant, const term_order *o) } template -gPoly::gPoly(const gSpace *p, const std::string &s, const term_order *o) - : Space(p), Order(o), Terms() +gPoly::gPoly(const VariableSpace *p, int var_no, int exp) + : Space(p) { - *this = s; // Operator = needs to be fixed + Terms.push_back(gMono((T)1, ExponentVector(p, var_no, exp))); } template -gPoly::gPoly(const gSpace *p, int var_no, int exp, const term_order *o) - : Space(p), Order(o), Terms() -{ - Terms.push_back(gMono((T)1, exp_vect(p, var_no, exp))); -} - -template -gPoly::gPoly(const gSpace *p, exp_vect exps, T coeff, const term_order *o) - : Space(p), Order(o), Terms() +gPoly::gPoly(const VariableSpace *p, ExponentVector exps, T coeff) + : Space(p) { Terms.push_back(gMono(coeff, exps)); } template -gPoly::gPoly(const gSpace *p, const gMono &mono, const term_order *o) - : Space(p), Order(o), Terms() +gPoly::gPoly(const VariableSpace *p, const gMono &mono) + : Space(p) { Terms.push_back(mono); } template -gPoly::gPoly(const gPoly &p) : Space(p.Space), Order(p.Order), Terms(p.Terms) +gPoly::gPoly(const gPoly &p) : Space(p.Space), Terms(p.Terms) { *this = p; } @@ -93,64 +85,6 @@ template gPoly &gPoly::operator=(const gPoly &p) return (*this); } -template gPoly &gPoly::operator=(const std::string &Hold) -{ - Gambit::List> nullTerms; - Terms = nullTerms; // get rid of old Terms - - charnum = 0; - int contflag = 1; - T nega = 1; - Gambit::Array PowArray(Space->Dmnsn()); - TheString = Hold + " +"; - - charc = TheString[charnum]; - - while (charnum <= TheString.length() && contflag) { - switch (charc) { - case '+': - case ' ': - charnum++; - charc = TheString[charnum]; - break; - case '-': - charnum++; - charc = TheString[charnum]; - nega = -nega; - break; - case 0: // Null termination of string - contflag = 0; - break; - default: - String_Term(nega); - nega = T(1); - break; - } - } - - Gambit::List> newTerms; - for (int j = 1; j <= Terms.size(); j++) { - int low = 0; - int high = newTerms.size() + 1; - while (low + 1 < high) { - int test = low + (high - low) / 2; - if (1 <= test && test <= newTerms.size()) { - // assert (Terms[j].ExpV() != newTerms[test].ExpV()); - } - if (Order->Less(Terms[j].ExpV(), newTerms[test].ExpV())) { - high = test; - } - else { - low = test; - } - } - newTerms.insert(std::next(newTerms.begin(), high - 1), Terms[j]); - } - Terms = newTerms; - - return (*this); -} - template gPoly gPoly::operator-() const { gPoly neg(*this); @@ -194,7 +128,7 @@ template void gPoly::operator+=(const gPoly &p) template void gPoly::operator+=(const T &val) { - *this += gPoly(Space, val, Order); + *this += gPoly(Space, val); } template gPoly gPoly::operator*(const gPoly &p) const @@ -246,181 +180,11 @@ template bool gPoly::operator==(const gPoly &p) const template bool gPoly::operator!=(const gPoly &p) const { return !(*this == p); } -//---------------------------------- -// Member Functions -//---------------------------------- - -template void gPoly::String_Term(T nega) -{ - Gambit::Array PowArray(Dmnsn()); - for (int a = 1; a <= Dmnsn(); a++) { - PowArray[a] = 0; - } - T val; - val = String_Coeff(nega); - - while (charc != '+' && charc != '-') { - if (charc == ' ') { - charnum++; - charc = TheString[charnum]; - } - else { - String_VarAndPow(PowArray); - } - } - - Terms.push_back(gMono(val, exp_vect(Space, PowArray))); -} - -template int gPoly::String_GetPow() -{ - - std::string Pow = ""; - while (charc == ' ') { - charnum++; - charc = TheString[charnum]; - } - - while (charc >= '0' && charc <= '9') { - Pow += charc; - charnum++; - charc = TheString[charnum]; - } - return (atoi(Pow.c_str())); -} - -template void gPoly::String_VarAndPow(Gambit::Array &PowArray) -{ - std::string VarName = ""; - int pow, varname; - while (charc != '^' && charc != ' ') { - VarName += charc; - charnum++; - charc = TheString[charnum]; - } - if (charc == '^') { - charnum++; - charc = TheString[charnum]; - pow = String_GetPow(); - } - else { - pow = 1; - } - for (varname = 1; varname <= Dmnsn() && VarName != (Space->VariableWithNumber(varname))->Name; - varname++) - ; - if (varname <= Dmnsn()) { - PowArray[varname] = pow; - } -} - -// bool function to check the string in &Hold - -template bool gPoly::Check_String(const std::string &Hold) -{ - unsigned int charnum = 0; - int boolflag = 0; - // state variables - int statenumber = 0; - int statevar = 0; - int statesign = 1; - // values of the state variables - enum number { nonumberbef, numberbef }; - enum var { novarbef, varbef }; - enum sign { nosignbef, signbef }; - std::string TheString = Hold; - char charc = TheString[charnum]; - // movement along the string with a switch case - while (charnum < TheString.length()) { - switch (charc) { - case '0': - case '1': - case '2': - case '3': - case '4': - case '5': - case '6': - case '7': - case '8': - case '9': - statenumber = numberbef; - statesign = nosignbef; - break; - case 'n': - if (statenumber == 0 && statesign == 0) { - boolflag++; - } - if (charnum == (TheString.length() - 1)) { - boolflag++; - } - statenumber = number(0); - statevar = varbef; - statesign = nosignbef; - break; - case '^': - if (statenumber == 0) { - boolflag++; - } - if (statevar == 0) { - boolflag++; - } - if (charnum == (TheString.length() - 1)) { - boolflag++; - } - statenumber = nonumberbef; - statevar = novarbef; - statesign = nosignbef; - break; - case '/': - case '.': - if (statenumber == 0) { - boolflag++; - } - if (charnum == (TheString.length() - 1)) { - boolflag++; - } - statenumber = nonumberbef; - statesign = nosignbef; - break; - case '+': - case '-': - if (statenumber == 0 && charnum != 0) { - boolflag++; - } - if (charnum == (TheString.length() - 1)) { - boolflag++; - } - statenumber = nonumberbef; - statevar = novarbef; - statesign = signbef; - break; - case ' ': - break; - default: - boolflag++; - break; - } - charnum++; - charc = TheString[charnum]; - } - // return values - if (boolflag == 0) { - return true; - } - else { - return false; - } -} - -template void gPoly::GetChar() { charc = TheString[charnum]; } - //---------------------------------- // Information //---------------------------------- -template const gSpace *gPoly::GetSpace() const { return (gSpace *)Space; } - -template const term_order *gPoly::GetOrder() const { return (term_order *)Order; } +template const VariableSpace *gPoly::GetSpace() const { return (VariableSpace *)Space; } template int gPoly::Dmnsn() const { return Space->Dmnsn(); } @@ -435,8 +199,6 @@ template int gPoly::DegreeOfVar(int var_no) const return max; } -template bool gPoly::IsZero() const { return Terms.empty(); } - template int gPoly::Degree() const { int max = 0; @@ -448,12 +210,12 @@ template int gPoly::Degree() const return max; } -template T gPoly::GetCoef(const Gambit::Array &Powers) const +template T gPoly::GetCoef(const Gambit::Vector &Powers) const { - return GetCoef(exp_vect(Space, Powers)); + return GetCoef(ExponentVector(Space, Powers)); } -template T gPoly::GetCoef(const exp_vect &Powers) const +template T gPoly::GetCoef(const ExponentVector &Powers) const { for (int j = 1; j <= Terms.size(); j++) { if (Terms[j].ExpV() == Powers) { @@ -493,7 +255,7 @@ template bool gPoly::IsMultiaffine() const return true; } -template T gPoly::Evaluate(const Gambit::Array &values) const +template T gPoly::Evaluate(const Gambit::Vector &values) const { T answer = 0; for (const auto &term : Terms) { @@ -502,57 +264,8 @@ template T gPoly::Evaluate(const Gambit::Array &values) const return answer; } -template Gambit::List gPoly::ExponentVectors() const -{ - Gambit::List result; - for (const auto &v : Terms) { - result.push_back(exp_vect(v.ExpV())); - } - return result; -} - template Gambit::List> gPoly::MonomialList() const { return Terms; } -template int gPoly::No_Monomials() const { return Terms.size(); } - -template int gPoly::UniqueActiveVariable() const -{ - Gambit::List ExpVecs = ExponentVectors(); - int activar = 0; - for (int i = 1; i <= ExpVecs.size(); i++) { - for (int j = 1; j <= Dmnsn(); j++) { - if (ExpVecs[i][j] > 0) { - if (activar > 0 && activar != j) { - return -1; // multivariate! - } - else { - activar = j; - } - } - } - } - return activar; -} - -template polynomial gPoly::UnivariateEquivalent(int activar) const -{ - // assert(UniqueActiveVariable() >= 0); - - Gambit::List coefs; - - if (!IsZero()) { - for (int h = 0; h <= DegreeOfVar(activar); h++) { - coefs.push_back((T)0); - } - - for (int i = 1; i <= Terms.size(); i++) { - coefs[Terms[i].ExpV()[activar] + 1] = Terms[i].Coef(); - } - } - - return polynomial(coefs); -} - //------------------------------------------------------------- // Private Versions of Arithmetic Operators //------------------------------------------------------------- @@ -582,11 +295,11 @@ Gambit::List> gPoly::Adder(const Gambit::List> &One, i++; } else { - if (Order->Less(One[i].ExpV(), Two[j].ExpV())) { + if (One[i].ExpV() < Two[j].ExpV()) { answer.push_back(One[i]); i++; } - else if (Order->Greater(One[i].ExpV(), Two[j].ExpV())) { + else if (One[i].ExpV() > Two[j].ExpV()) { answer.push_back(Two[j]); j++; } @@ -625,10 +338,10 @@ Gambit::List> gPoly::Mult(const Gambit::List> &One, else { int bot = 1; int top = answer.size(); - if (Order->Less(answer[top].ExpV(), next.ExpV())) { + if (answer[top].ExpV() < next.ExpV()) { answer.push_back(next); } - else if (Order->Greater(answer[bot].ExpV(), next.ExpV())) { + else if (answer[bot].ExpV() > next.ExpV()) { answer.push_front(next); } else { @@ -644,10 +357,10 @@ Gambit::List> gPoly::Mult(const Gambit::List> &One, if (answer[test].ExpV() == next.ExpV()) { bot = top = test; } - else if (Order->Less(answer[test].ExpV(), next.ExpV())) { + else if (answer[test].ExpV() < next.ExpV()) { bot = test; } - else { // (Order->Greater(answer[test].ExpV(),next.ExpV())) + else { // (answer[test].ExpV() > next.ExpV()) top = test; } } @@ -667,13 +380,11 @@ Gambit::List> gPoly::Mult(const Gambit::List> &One, template gPoly gPoly::DivideByPolynomial(const gPoly &den) const { - gPoly zero(Space, (T)0, Order); + gPoly zero(Space, (T)0); if (den == zero) { throw Gambit::ZeroDivideException(); } - // assert(*this == zero || den.Degree() <= Degree()); - // assumes exact divisibility! gPoly result = zero; @@ -695,8 +406,7 @@ template gPoly gPoly::DivideByPolynomial(const gPoly &den) co while (remainder != zero) { gPoly quot = remainder.LeadingCoefficient(last) / den.LeadingCoefficient(last); - gPoly power_of_last(Space, last, remainder.DegreeOfVar(last) - den.DegreeOfVar(last), - Order); + gPoly power_of_last(Space, last, remainder.DegreeOfVar(last) - den.DegreeOfVar(last)); result += quot * power_of_last; remainder -= quot * power_of_last * den; } @@ -704,39 +414,13 @@ template gPoly gPoly::DivideByPolynomial(const gPoly &den) co return result; } -template gPoly gPoly::EvaluateOneVar(int varnumber, T val) const -{ - gPoly answer(Space, (T)0, Order); - - for (int i = 1; i <= Terms.size(); i++) { - answer += - gPoly(Space, Terms[i].ExpV().AfterZeroingOutExpOfVariable(varnumber), - ((T)Terms[i].Coef()) * ((T)pow(val, (double)Terms[i].ExpV()[varnumber])), Order); - } - return answer; -} - -template -exp_vect gPoly::OrderMaxMonomialDivisibleBy(const term_order &order, const exp_vect & /*expv*/) -{ - // gout << "You have just tested OrderMaxMonomialDivisibleBy.\n"; - - exp_vect answer(Space); // constructs [0,..,0] - for (int i = 1; i <= Terms.size(); i++) { - if (order.Less(answer, Terms[i].ExpV()) && answer < Terms[i].ExpV()) { - answer = Terms[i].ExpV(); - } - } - return answer; -} - template gPoly gPoly::PartialDerivative(int varnumber) const { gPoly newPoly(*this); for (int i = 1; i <= newPoly.Terms.size(); i++) { newPoly.Terms[i] = gMono(newPoly.Terms[i].Coef() * (T)newPoly.Terms[i].ExpV()[varnumber], - newPoly.Terms[i].ExpV().AfterDecrementingExpOfVariable(varnumber)); + newPoly.Terms[i].ExpV().WithZeroExponent(varnumber)); } int j = 1; @@ -762,117 +446,22 @@ template gPoly gPoly::LeadingCoefficient(int varnumber) const for (int j = 1; j <= Terms.size(); j++) { if (Terms[j].ExpV()[varnumber] == degree) { newPoly.Terms.push_back( - gMono(Terms[j].Coef(), Terms[j].ExpV().AfterZeroingOutExpOfVariable(varnumber))); + gMono(Terms[j].Coef(), Terms[j].ExpV().WithZeroExponent(varnumber))); } } return (newPoly); } -//-------------------- -// Term Order Concepts -//-------------------- - -template exp_vect gPoly::LeadingPowerProduct(const term_order &order) const -{ - // assert (Terms.Length() > 0); - - if (*Order == order) { // worth a try ... - return Terms.back().ExpV(); - } - else { - int max = 1; - for (int j = 2; j <= Terms.size(); j++) { - if (order.Less(Terms[max].ExpV(), Terms[j].ExpV())) { - max = j; - } - } - return Terms[max].ExpV(); - } -} - -template T gPoly::LeadingCoefficient(const term_order &order) const -{ - if (*Order == order) { // worth a try ... - return Terms.back().Coef(); - } - else { - int max = 1; - for (int j = 2; j <= Terms.size(); j++) { - if (order.Less(Terms[max].ExpV(), Terms[j].ExpV())) { - max = j; - } - } - return Terms[max].Coef(); - } -} - -template gPoly gPoly::LeadingTerm(const term_order &order) const -{ - if (*Order == order) { // worth a try ... - return gPoly(Space, Terms.back(), Order); - } - else { - int max = 1; - for (int j = 2; j <= Terms.size(); j++) { - if (order.Less(Terms[max].ExpV(), Terms[j].ExpV())) { - max = j; - } - } - return gPoly(Space, Terms[max], Order); - } -} - -template void gPoly::ToMonic(const term_order &order) -{ - *this = *this / LeadingCoefficient(order); -} - -template -void gPoly::ReduceByDivisionAtExpV(const term_order &order, const gPoly &divisor, - const exp_vect &expv) -{ - // assert(expv >= divisor.LeadingPowerProduct(order)); - // assert(divisor.LeadingCoefficient(order) != (T)0); - - gPoly factor(Space, expv - divisor.LeadingPowerProduct(order), (T)1, Order); - - *this -= (GetCoef(expv) / divisor.LeadingCoefficient(order)) * factor * divisor; -} - -template -void gPoly::ReduceByRepeatedDivision(const term_order &order, const gPoly &divisor) -{ - exp_vect zero_exp_vec(Space); - - exp_vect exps = OrderMaxMonomialDivisibleBy(order, divisor.LeadingPowerProduct(order)); - - while (exps != zero_exp_vec) { - ReduceByDivisionAtExpV(order, divisor, exps); - exps = OrderMaxMonomialDivisibleBy(order, divisor.LeadingPowerProduct(order)); - } -} - -template -gPoly gPoly::S_Polynomial(const term_order &order, const gPoly &arg2) const -{ - exp_vect exp_lcm = (LeadingPowerProduct(order)).LCM(arg2.LeadingPowerProduct(order)); - gPoly lcm = gPoly(Space, exp_lcm, (T)1, Order); - gPoly L1 = lcm / LeadingTerm(order); - gPoly L2 = lcm / arg2.LeadingTerm(order); - - return L1 * (*this) - L2 * arg2; -} - template gPoly gPoly::TranslateOfMono(const gMono &m, const Gambit::Vector &new_origin) const { - gPoly answer(GetSpace(), (T)1, GetOrder()); + gPoly answer(GetSpace(), (T)1); for (int i = 1; i <= Dmnsn(); i++) { if (m.ExpV()[i] > 0) { - gPoly lt(GetSpace(), i, 1, GetOrder()); - lt += gPoly(GetSpace(), new_origin[i], GetOrder()); + gPoly lt(GetSpace(), i, 1); + lt += gPoly(GetSpace(), new_origin[i]); for (int j = 1; j <= m.ExpV()[i]; j++) { answer *= lt; } @@ -886,7 +475,7 @@ gPoly gPoly::TranslateOfMono(const gMono &m, const Gambit::Vector &n template gPoly gPoly::TranslateOfPoly(const Gambit::Vector &new_origin) const { - gPoly answer(GetSpace(), GetOrder()); + gPoly answer(GetSpace()); for (int i = 1; i <= this->MonomialList().size(); i++) { answer += TranslateOfMono(this->MonomialList()[i], new_origin); } @@ -898,14 +487,14 @@ gPoly gPoly::MonoInNewCoordinates(const gMono &m, const Gambit::SquareM { // assert(M.NumRows() == Dmnsn()); - gPoly answer(GetSpace(), (T)1, GetOrder()); + gPoly answer(GetSpace(), (T)1); for (int i = 1; i <= Dmnsn(); i++) { if (m.ExpV()[i] > 0) { - gPoly linearform(GetSpace(), (T)0, GetOrder()); + gPoly linearform(GetSpace(), (T)0); for (int j = 1; j <= Dmnsn(); j++) { - exp_vect exps(GetSpace(), j, 1); - linearform += gPoly(GetSpace(), exps, M(i, j), GetOrder()); + ExponentVector exps(GetSpace(), j, 1); + linearform += gPoly(GetSpace(), exps, M(i, j)); } for (int k = 1; k <= m.ExpV()[i]; k++) { answer *= linearform; @@ -920,7 +509,7 @@ gPoly gPoly::MonoInNewCoordinates(const gMono &m, const Gambit::SquareM template gPoly gPoly::PolyInNewCoordinates(const Gambit::SquareMatrix &M) const { - gPoly answer(GetSpace(), GetOrder()); + gPoly answer(GetSpace()); for (int i = 1; i <= MonomialList().size(); i++) { answer += MonoInNewCoordinates(MonomialList()[i], M); } @@ -961,74 +550,18 @@ template gPoly operator+(const T &val, const gPoly &poly) template gPoly operator+(const gPoly &poly, const T &val) { return val + poly; } -template void gPoly::Output(std::string &t) const -{ - std::string s; - if (Terms.empty()) { - s += "0"; - } - else { - for (int j = 1; j <= Terms.size(); j++) { - if (Terms[j].Coef() < (T)0) { - s += "-"; - if (j > 1) { - s += " "; - } - } - else if (j > 1) { - s += "+ "; - } - - if ((Terms[j].Coef() != (T)1 && -Terms[j].Coef() != (T)1) || Terms[j].IsConstant()) { - if (Terms[j].Coef() < (T)0) { - s += Gambit::lexical_cast(-Terms[j].Coef()); - } - else { - s += Gambit::lexical_cast(Terms[j].Coef()); - } - - for (int k = 1; k <= Space->Dmnsn(); k++) { - int exp = Terms[j].ExpV()[k]; - if (exp > 0) { - s += " "; - s += (*Space)[k]->Name; - if (exp != 1) { - s += '^'; - s += Gambit::lexical_cast(exp); - } - } - } - - if (j < Terms.size()) { - s += " "; - } - } - } - } - if (s == "") { - s = " 0"; - } - - t += s; -} - -template std::string &operator<<(std::string &p_text, const gPoly &p_poly) -{ - p_poly.Output(p_text); - return p_text; -} //---------------------------------- // Conversion //---------------------------------- -template gPoly TogDouble(const gPoly &given) +template gPoly ToDouble(const gPoly &given) { - gPoly answer(given.GetSpace(), given.GetOrder()); + gPoly answer(given.GetSpace()); Gambit::List> list = given.MonomialList(); for (int i = 1; i <= list.size(); i++) { auto nextcoef = (double)list[i].Coef(); - gPoly next(given.GetSpace(), list[i].ExpV(), nextcoef, given.GetOrder()); + gPoly next(given.GetSpace(), list[i].ExpV(), nextcoef); answer += next; } @@ -1044,13 +577,13 @@ template gPoly NormalizationOfPoly(const gPoly &given) } if (maxcoeff < 0.000000001) { - return TogDouble(given); + return ToDouble(given); } - gPoly answer(given.GetSpace(), given.GetOrder()); + gPoly answer(given.GetSpace()); for (int i = 1; i <= list.size(); i++) { double nextcoef = (double)list[i].Coef() / maxcoeff; - gPoly next(given.GetSpace(), list[i].ExpV(), nextcoef, given.GetOrder()); + gPoly next(given.GetSpace(), list[i].ExpV(), nextcoef); answer += next; } diff --git a/src/solvers/enumpoly/gpolylst.h b/src/solvers/enumpoly/gpolylst.h index 6569c7f78..87a98e20f 100644 --- a/src/solvers/enumpoly/gpolylst.h +++ b/src/solvers/enumpoly/gpolylst.h @@ -29,15 +29,14 @@ template class gPolyList { private: - const gSpace *Space; - const term_order *Order; + const VariableSpace *Space; Gambit::List *> List; public: - gPolyList(const gSpace *sp, const term_order *to) : Space(sp), Order(to) {} + gPolyList(const VariableSpace *sp) : Space(sp) {} - gPolyList(const gSpace *, const term_order *, const Gambit::List *> &); - gPolyList(const gSpace *, const term_order *, const Gambit::List> &); + gPolyList(const VariableSpace *, const Gambit::List *> &); + gPolyList(const VariableSpace *, const Gambit::List> &); gPolyList(const gPolyList &); ~gPolyList(); // Deletes all pointees @@ -60,8 +59,7 @@ template class gPolyList { gPolyList SystemInNewCoordinates(const Gambit::SquareMatrix &) const; // Information - const gSpace *AmbientSpace() const { return Space; } - const term_order *TermOrder() const { return Order; } + const VariableSpace *AmbientSpace() const { return Space; } int Length() const { return List.size(); } int Dmnsn() const { return Space->Dmnsn(); } bool IsMultiaffine() const; diff --git a/src/solvers/enumpoly/gpolylst.imp b/src/solvers/enumpoly/gpolylst.imp index 06aef6527..e63d9562c 100644 --- a/src/solvers/enumpoly/gpolylst.imp +++ b/src/solvers/enumpoly/gpolylst.imp @@ -44,9 +44,8 @@ Gambit::List InteriorSegment(const Gambit::List &p_list, int first, int la template -gPolyList::gPolyList(const gSpace *sp, const term_order *to, - const Gambit::List *> &plist) - : Space(sp), Order(to), List() +gPolyList::gPolyList(const VariableSpace *sp, const Gambit::List *> &plist) + : Space(sp) { for (int ii = 1; ii <= plist.size(); ii++) { auto *temp = new gPoly(*plist[ii]); @@ -55,8 +54,8 @@ gPolyList::gPolyList(const gSpace *sp, const term_order *to, } template -gPolyList::gPolyList(const gSpace *sp, const term_order *to, const Gambit::List> &list) - : Space(sp), Order(to), List() +gPolyList::gPolyList(const VariableSpace *sp, const Gambit::List> &list) + : Space(sp) { for (int ii = 1; ii <= list.size(); ii++) { auto *temp = new gPoly(list[ii]); @@ -65,7 +64,7 @@ gPolyList::gPolyList(const gSpace *sp, const term_order *to, const Gambit::Li } template -gPolyList::gPolyList(const gPolyList &lst) : Space(lst.Space), Order(lst.Order), List() +gPolyList::gPolyList(const gPolyList &lst) : Space(lst.Space) { for (int ii = 1; ii <= lst.List.size(); ii++) { auto *temp = new gPoly(*(lst.List[ii])); @@ -102,7 +101,7 @@ template gPolyList &gPolyList::operator=(const gPolyList &rhs template bool gPolyList::operator==(const gPolyList &rhs) const { - if (Space != rhs.Space || Order != rhs.Order) { + if (Space != rhs.Space) { return false; } if (List.size() != rhs.List.size()) { @@ -142,7 +141,7 @@ gPolyList gPolyList::TranslateOfSystem(const Gambit::Vector &new_origin for (int i = 1; i <= Length(); i++) { new_polys.push_back((*this)[i].TranslateOfPoly(new_origin)); } - return gPolyList(AmbientSpace(), TermOrder(), new_polys); + return gPolyList(AmbientSpace(), new_polys); } template @@ -152,7 +151,7 @@ gPolyList gPolyList::SystemInNewCoordinates(const Gambit::SquareMatrix for (int i = 1; i <= Length(); i++) { new_polys.push_back((*this)[i].PolyInNewCoordinates(M)); } - return gPolyList(AmbientSpace(), TermOrder(), new_polys); + return gPolyList(AmbientSpace(), new_polys); } //---------------------------------- diff --git a/src/solvers/enumpoly/ineqsolv.h b/src/solvers/enumpoly/ineqsolv.h index 5ce2f5bbc..328e2ac9b 100644 --- a/src/solvers/enumpoly/ineqsolv.h +++ b/src/solvers/enumpoly/ineqsolv.h @@ -75,8 +75,7 @@ template class IneqSolv { bool operator!=(const IneqSolv &) const; // Information - const gSpace *AmbientSpace() const { return System.AmbientSpace(); } - const term_order *TermOrder() const { return System.TermOrder(); } + const VariableSpace *AmbientSpace() const { return System.AmbientSpace(); } int Dmnsn() const { return System.Dmnsn(); } const gPolyList &UnderlyingEquations() const { return System; } T ErrorTolerance() const { return Epsilon; } diff --git a/src/solvers/enumpoly/monomial.cc b/src/solvers/enumpoly/monomial.cc deleted file mode 100644 index 38bae6b99..000000000 --- a/src/solvers/enumpoly/monomial.cc +++ /dev/null @@ -1,27 +0,0 @@ -// -// This file is part of Gambit -// Copyright (c) 1994-2024, The Gambit Project (http://www.gambit-project.org) -// -// FILE: src/tools/enumpoly/monomial.cc -// Instantiation of monomial classes -// -// 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 "monomial.imp" - -template class gMono; -template class gMono; -// template class gMono; diff --git a/src/solvers/enumpoly/monomial.h b/src/solvers/enumpoly/monomial.h deleted file mode 100644 index d40dd9ede..000000000 --- a/src/solvers/enumpoly/monomial.h +++ /dev/null @@ -1,66 +0,0 @@ -// -// This file is part of Gambit -// Copyright (c) 1994-2024, The Gambit Project (http://www.gambit-project.org) -// -// FILE: src/tools/enumpoly/monomial.h -// Declaration of monomial class -// -// 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 "prepoly.h" - -// This file provides the template class -// -// gMono -// -// whose objects are monomials in several variables -// with coefficients of class T and nonnegative exponents. -// This role of this class is to support the class gPoly. - -template class gMono { -private: - T coef; - exp_vect exps; - -public: - // constructors - gMono(const gSpace *, const T &); - gMono(const T &, const exp_vect &); - gMono(const gMono &); - ~gMono(); - - // operators - gMono &operator=(const gMono &); - - bool operator==(const gMono &) const; - bool operator!=(const gMono &) const; - gMono operator*(const gMono &) const; - gMono operator/(const gMono &) const; - gMono operator+(const gMono &) const; // assert exps == - gMono &operator+=(const gMono &); // assert exps == - gMono &operator*=(const T &); - gMono operator-() const; - - // information - const T &Coef() const; - int Dmnsn() const; - int TotalDegree() const; - bool IsConstant() const; - bool IsMultiaffine() const; - const exp_vect &ExpV() const; - T Evaluate(const Gambit::Array &) const; - T Evaluate(const Gambit::Vector &) const; -}; diff --git a/src/solvers/enumpoly/monomial.imp b/src/solvers/enumpoly/monomial.imp deleted file mode 100644 index 0214402e4..000000000 --- a/src/solvers/enumpoly/monomial.imp +++ /dev/null @@ -1,134 +0,0 @@ -// -// This file is part of Gambit -// Copyright (c) 1994-2024, The Gambit Project (http://www.gambit-project.org) -// -// FILE: src/tools/enumpoly/monomial.imp -// Implementation of monomial classes -// -// 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 "monomial.h" - -//-------------------------------------------------------------------------- -// gMono -- constructors and destructor -//-------------------------------------------------------------------------- - -template gMono::gMono(const gSpace *p, const T &x) : coef(x), exps(p) {} - -template gMono::gMono(const T &x, const exp_vect &e) : coef(x), exps(e) -{ - if (x == (T)0) { - exps.ToZero(); - } -} - -template gMono::gMono(const gMono &y) : coef(y.coef), exps(y.exps) {} - -template gMono::~gMono() = default; - -//-------------------------------------------------------------------------- -// gMono -- operators -//-------------------------------------------------------------------------- - -template gMono &gMono::operator=(const gMono &y) -{ - if (this != &y) { - coef = y.coef; - exps = y.exps; - } - return *this; -} - -template bool gMono::operator==(const gMono &y) const -{ - return (coef == y.coef && exps == y.exps); -} - -template bool gMono::operator!=(const gMono &y) const { return !(*this == y); } - -template gMono gMono::operator*(const gMono &y) const -{ - return gMono(coef * y.coef, exps + y.exps); -} - -template gMono gMono::operator/(const gMono &y) const -{ - // assert ( y.coef != (T)0); - return gMono(coef / y.coef, exps - y.exps); -} - -template gMono gMono::operator+(const gMono &y) const -{ - // assert (exps == y.exps); - return gMono(coef + y.coef, exps); -} - -template gMono &gMono::operator+=(const gMono &y) -{ - // assert (exps == y.exps); - coef += y.coef; - return *this; -} - -template gMono &gMono::operator*=(const T &val) -{ - coef *= val; - return *this; -} - -template gMono gMono::operator-() const { return gMono(-coef, exps); } - -//-------------------------------------------------------------------------- -// gMono -- information -//-------------------------------------------------------------------------- - -template const T &gMono::Coef() const { return coef; } - -template int gMono::Dmnsn() const { return exps.Dmnsn(); } - -template int gMono::TotalDegree() const { return exps.TotalDegree(); } - -template const exp_vect &gMono::ExpV() const { return exps; } - -template bool gMono::IsConstant() const { return exps.IsConstant(); } - -template bool gMono::IsMultiaffine() const { return exps.IsMultiaffine(); } - -template T gMono::Evaluate(const Gambit::Array &vals) const -{ - T answer = Coef(); - - for (int i = 1; i <= Dmnsn(); i++) { - for (int j = 1; j <= exps[i]; j++) { - answer *= vals[i]; - } - } - - return answer; -} - -template T gMono::Evaluate(const Gambit::Vector &vals) const -{ - T answer = Coef(); - - for (int i = 1; i <= Dmnsn(); i++) { - for (int j = 1; j <= exps[i]; j++) { - answer *= vals[i]; - } - } - - return answer; -} diff --git a/src/solvers/enumpoly/nfgpoly.cc b/src/solvers/enumpoly/nfgpoly.cc index 8ff23a99a..3521f9dfd 100644 --- a/src/solvers/enumpoly/nfgpoly.cc +++ b/src/solvers/enumpoly/nfgpoly.cc @@ -37,18 +37,17 @@ namespace { // The polynomial representation of each strategy probability, substituting in // the sum-to-one equation for the probability of the last strategy for each player -std::map> BuildStrategyVariables(const gSpace &space, - const term_order &lex, +std::map> BuildStrategyVariables(const VariableSpace &space, const StrategySupportProfile &support) { int index = 1; std::map> strategy_poly; for (auto player : support.GetGame()->GetPlayers()) { auto strategies = support.GetStrategies(player); - gPoly residual(&space, 1, &lex); + gPoly residual(&space, 1); for (auto strategy : strategies) { if (strategy != strategies.back()) { - strategy_poly.try_emplace(strategy, &space, index, 1, &lex); + strategy_poly.try_emplace(strategy, &space, index, 1); residual -= strategy_poly.at(strategy); index++; } @@ -60,15 +59,15 @@ std::map> BuildStrategyVariables(const gSpace &space return strategy_poly; } -gPoly IndifferenceEquation(const gSpace &space, const term_order &lex, +gPoly IndifferenceEquation(const VariableSpace &space, const StrategySupportProfile &support, const std::map> &strategy_poly, const GameStrategy &s1, const GameStrategy &s2) { - gPoly equation(&space, &lex); + gPoly equation(&space); for (auto iter : StrategyContingencies(support, {s1})) { - gPoly term(&space, 1, &lex); + gPoly term(&space, 1); for (auto player : support.GetGame()->GetPlayers()) { if (player != s1->GetPlayer()) { term *= strategy_poly.at(iter->GetStrategy(player)); @@ -80,17 +79,17 @@ gPoly IndifferenceEquation(const gSpace &space, const term_order &lex, return equation; } -gPolyList ConstructEquations(const gSpace &space, const term_order &lex, +gPolyList ConstructEquations(const VariableSpace &space, const StrategySupportProfile &support, const std::map> &strategy_poly) { - gPolyList equations(&space, &lex); + gPolyList equations(&space); // Indifference equations between pairs of strategies for each player for (auto player : support.GetPlayers()) { auto strategies = support.GetStrategies(player); for (auto s1 = strategies.begin(), s2 = std::next(strategies.begin()); s2 != strategies.end(); ++s1, ++s2) { - equations += IndifferenceEquation(space, lex, support, strategy_poly, *s1, *s2); + equations += IndifferenceEquation(space, support, strategy_poly, *s1, *s2); } } // Inequalities for last probability for each player @@ -109,11 +108,10 @@ std::list> EnumPolyStrategySupportSolve(const StrategySupportProfile &support, bool &is_singular, int p_stopAfter) { - gSpace Space(support.MixedProfileLength() - support.GetGame()->NumPlayers()); - term_order Lex(&Space, lex); + VariableSpace Space(support.MixedProfileLength() - support.GetGame()->NumPlayers()); - auto strategy_poly = BuildStrategyVariables(Space, Lex, support); - gPolyList equations = ConstructEquations(Space, Lex, support, strategy_poly); + auto strategy_poly = BuildStrategyVariables(Space, support); + gPolyList equations = ConstructEquations(Space, support, strategy_poly); Vector bottoms(Space.Dmnsn()), tops(Space.Dmnsn()); bottoms = 0; diff --git a/src/solvers/enumpoly/poly.cc b/src/solvers/enumpoly/poly.cc deleted file mode 100644 index 5e1fb0840..000000000 --- a/src/solvers/enumpoly/poly.cc +++ /dev/null @@ -1,26 +0,0 @@ -// -// This file is part of Gambit -// Copyright (c) 1994-2024, The Gambit Project (http://www.gambit-project.org) -// -// FILE: src/tools/enumpoly/poly.cc -// Instantiation of polynomial classes -// -// 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 "poly.imp" - -template class polynomial; -template class polynomial; diff --git a/src/solvers/enumpoly/poly.h b/src/solvers/enumpoly/poly.h deleted file mode 100644 index d6add239f..000000000 --- a/src/solvers/enumpoly/poly.h +++ /dev/null @@ -1,150 +0,0 @@ -// -// This file is part of Gambit -// Copyright (c) 1994-2024, The Gambit Project (http://www.gambit-project.org) -// -// FILE: src/tools/enumpoly/poly.h -// Declaration of polynomial classes -// -// 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 "gambit.h" -#include "interval.h" -#include "gcomplex.h" - -/* This file supplies the template class - - polynomial - -These are univariate polynomials with coefficients of class T. -Polynomials are implemented as Gambit::List's of coefficients. There is no -attempt to maintain sparseness. - -*/ - -template class polynomial { - -private: - Gambit::List coeflist; - -public: - // constructors and destructor - explicit polynomial(int = -1); - polynomial(const polynomial &); - explicit polynomial(const Gambit::List &); - explicit polynomial(const Gambit::Vector &); - polynomial(const T &, const int &); - ~polynomial() = default; - - // unary operators - polynomial operator-() const; - polynomial Derivative() const; - - // binary operators - polynomial &operator=(const polynomial &y); - bool operator==(const polynomial &y) const; - bool operator!=(const polynomial &y) const; - const T &operator[](int index) const; - polynomial operator+(const polynomial &y) const; - polynomial operator-(const polynomial &y) const; - polynomial operator*(const polynomial &y) const; - polynomial operator/(const polynomial &y) const; - polynomial &operator+=(const polynomial &y); - polynomial &operator-=(const polynomial &y); - polynomial &operator*=(const polynomial &y); - polynomial &operator/=(const polynomial &y); - polynomial operator%(const polynomial &y) const; - - // manipulation - void ToMonic(); - // polynomial Togdouble() const; - - polynomial TogDouble() const; - - // information - bool IsZero() const; - T EvaluationAt(const T &arg) const; - int Degree() const; - T LeadingCoefficient() const; - Gambit::List CoefficientList() const; - polynomial GcdWith(const polynomial &) const; - bool IsQuadratfrei() const; - bool CannotHaveRootsIn(const gInterval &) const; - Gambit::List> RootSubintervals(const gInterval &) const; - gInterval NeighborhoodOfRoot(const gInterval &, T &) const; - Gambit::List> PreciseRootIntervals(const gInterval &, T &) const; - Gambit::List PreciseRoots(const gInterval &, T &) const; -}; - -/* REMARKS - - The function cannot_have_roots_in is based on the principle that if - -f = a_0 + a_1x + ... + a_dx^d - -with a_0 > 0, then - -abs(f(t)) >= a_0 - max{abs(a_1),...,abs(a_d)}*(abs(t) + ... + abs(t)^d) - -and the RHS will be positive whenever - -//WRONG! abs(t) < a_0/(a_0 + max{abs(a_1),...,abs(a_d)}). - -*/ - -class complexpoly { - -private: - Gambit::List coeflist; - -public: - // constructors and destructor - explicit complexpoly(int = -1); - complexpoly(const complexpoly &); - explicit complexpoly(const Gambit::List &); - complexpoly(const gComplex &, const int &); - ~complexpoly(); - - // unary operators - complexpoly operator-() const; - complexpoly Derivative() const; - - // binary operators - complexpoly &operator=(const complexpoly &y); - bool operator==(const complexpoly &y) const; - bool operator!=(const complexpoly &y) const; - const gComplex &operator[](int index) const; - complexpoly operator+(const complexpoly &y) const; - complexpoly operator-(const complexpoly &y) const; - complexpoly operator*(const complexpoly &y) const; - complexpoly operator/(const complexpoly &y) const; - complexpoly &operator+=(const complexpoly &y); - complexpoly &operator-=(const complexpoly &y); - complexpoly &operator*=(const complexpoly &y); - complexpoly &operator/=(const complexpoly &y); - complexpoly operator%(const complexpoly &y) const; - - // manipulation - void ToMonic(); - - // information - bool IsZero() const; - gComplex EvaluationAt(const gComplex &arg) const; - int Degree() const; - gComplex LeadingCoefficient() const; - complexpoly GcdWith(const complexpoly &) const; - bool IsQuadratfrei() const; - Gambit::List Roots() const; -}; diff --git a/src/solvers/enumpoly/poly.imp b/src/solvers/enumpoly/poly.imp deleted file mode 100644 index 3629965f6..000000000 --- a/src/solvers/enumpoly/poly.imp +++ /dev/null @@ -1,862 +0,0 @@ -// -// This file is part of Gambit -// Copyright (c) 1994-2024, The Gambit Project (http://www.gambit-project.org) -// -// FILE: src/tools/enumpoly/poly.imp -// Implementation of polynomial class -// -// 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" - -//-------------------------------------------------------------------------- -// class: polynomial -//-------------------------------------------------------------------------- - -//-------------------------------------------------------------------------- -// constructors and a destructor -//-------------------------------------------------------------------------- - -template polynomial::polynomial(const polynomial &x) : coeflist(x.coeflist) {} - -template polynomial::polynomial(const T &coeff, const int °) -{ - if (coeff != (T)0) { - for (int i = 0; i < deg; i++) { - coeflist.push_back((T)0); - } - coeflist.push_back(coeff); - } -} - -template -polynomial::polynomial(const Gambit::List &coefficientlist) : coeflist(coefficientlist) -{ -} - -template -polynomial::polynomial(const Gambit::Vector &coefficientvector) : coeflist() -{ - for (int i = 1; i <= coefficientvector.Length(); i++) { - coeflist.push_back(coefficientvector[i]); - } -} - -template polynomial::polynomial(int deg) : coeflist() -{ - if (deg >= 0) { - // gout << "Error is polynomial int constructor.\n"; - exit(1); - } -} - -//-------------------------------------------------------------------------- -// operators -//-------------------------------------------------------------------------- - -template polynomial &polynomial::operator=(const polynomial &y) -{ - if (this != &y) { - coeflist = y.coeflist; - } - - return *this; -} - -template bool polynomial::operator==(const polynomial &y) const -{ - if (Degree() != y.Degree()) { - return false; - } - else { - for (int i = 0; i <= Degree(); i++) { - if (coeflist[i + 1] != y.coeflist[i + 1]) { - return false; - } - } - } - return true; -} - -template bool polynomial::operator!=(const polynomial &y) const -{ - return !(*this == y); -} - -template const T &polynomial::operator[](const int index) const -{ - return coeflist[index + 1]; -} - -template polynomial polynomial::operator+(const polynomial &y) const -{ - if (Degree() < 0) { - return polynomial(y); - } - else if (y.Degree() < 0) { - return polynomial(*this); - } - - int max_degree; - - if (Degree() > y.Degree()) { - max_degree = Degree(); - } - else { - max_degree = y.Degree(); - } - - polynomial sum; - for (int i = 0; i <= max_degree; i++) { - sum.coeflist.push_back((T)0); - if (i <= Degree()) { - sum.coeflist[i + 1] += coeflist[i + 1]; - } - if (i <= y.Degree()) { - sum.coeflist[i + 1] += y.coeflist[i + 1]; - } - } - - while (!sum.coeflist.empty() && sum.coeflist.back() == (T)0) { - sum.coeflist.pop_back(); - } - - return sum; -} - -template polynomial polynomial::operator-(const polynomial &y) const -{ - return polynomial(*this + (-y)); -} - -template polynomial polynomial::operator*(const polynomial &y) const -{ - if (Degree() == -1) { - return polynomial(*this); - } - else if (y.Degree() == -1) { - return polynomial(y); - } - - int tot_degree = Degree() + y.Degree(); - - polynomial product; - for (int t = 0; t <= tot_degree; t++) { - product.coeflist.push_back((T)0); - } - for (int i = 0; i <= Degree(); i++) { - for (int j = 0; j <= y.Degree(); j++) { - product.coeflist[i + j + 1] += (*this)[i] * y[j]; - } - } - - return product; -} - -template polynomial polynomial::operator/(const polynomial &q) const -{ - // assert(q.Degree() >= 0); - - polynomial ans; - polynomial r = *this; - while (r.Degree() >= q.Degree()) { - polynomial x(r.LeadingCoefficient() / q.LeadingCoefficient(), r.Degree() - q.Degree()); - ans += x; - r -= q * x; - } - return polynomial(ans); -} - -template polynomial &polynomial::operator+=(const polynomial &y) -{ - return ((*this) = (*this) + y); -} - -template polynomial &polynomial::operator-=(const polynomial &y) -{ - return ((*this) = (*this) - y); -} - -template polynomial &polynomial::operator*=(const polynomial &y) -{ - return ((*this) = (*this) * y); -} - -template polynomial &polynomial::operator/=(const polynomial &y) -{ - return ((*this) = (*this) / y); -} - -template polynomial polynomial::operator%(const polynomial &q) const -{ - // assert (q.Degree() != -1); - - polynomial ans; - polynomial r = *this; - - while (r.Degree() >= q.Degree()) { - polynomial x(r.LeadingCoefficient() / q.LeadingCoefficient(), r.Degree() - q.Degree()); - ans += x; - r -= q * x; - } - return polynomial(r); -} - -template polynomial polynomial::operator-() const -{ - polynomial negation; - for (int i = 0; i <= Degree(); i++) { - negation.coeflist.push_back(-coeflist[i + 1]); - } - - return negation; -} - -template polynomial polynomial::Derivative() const -{ - if (Degree() <= 0) { - return polynomial(-1); - } - - polynomial derivative; - for (int i = 1; i <= Degree(); i++) { - derivative.coeflist.push_back((T)i * coeflist[i + 1]); - } - - return derivative; -} - -//-------------------------------------------------------------------------- -// manipulation -//-------------------------------------------------------------------------- - -template void polynomial::ToMonic() -{ - // assert (!IsZero()); - - T lc = LeadingCoefficient(); - for (int i = 1; i <= coeflist.size(); i++) { - coeflist[i] /= lc; - } -} - -template polynomial polynomial::TogDouble() const -{ - Gambit::List newcoefs; - for (const auto &coef : coeflist) { - newcoefs.push_back((double)coef); - } - return polynomial(newcoefs); -} - -//-------------------------------------------------------------------------- -// information -//-------------------------------------------------------------------------- - -template bool polynomial::IsZero() const { return coeflist.empty(); } - -template T polynomial::EvaluationAt(const T &arg) const -{ - T answer; - - if (IsZero()) { - answer = (T)0; - } - else { - answer = coeflist[Degree() + 1]; - for (int i = Degree(); i >= 1; i--) { - answer *= arg; - answer += coeflist[i]; - } - } - - return answer; -} - -template int polynomial::Degree() const { return coeflist.size() - 1; } - -template T polynomial::LeadingCoefficient() const -{ - if (Degree() < 0) { - return (T)0; - } - else { - return coeflist[Degree() + 1]; - } -} - -template Gambit::List polynomial::CoefficientList() const { return coeflist; } - -template polynomial polynomial::GcdWith(const polynomial &that) const -{ - // assert( !this->IsZero() && !that.IsZero() ); - - polynomial numerator(*this); - numerator.ToMonic(); - polynomial denominator(that); - denominator.ToMonic(); - polynomial remainder(numerator % denominator); - - while (!remainder.IsZero()) { - remainder.ToMonic(); - numerator = denominator; - denominator = remainder; - remainder = numerator % denominator; - } - - return denominator; -} - -template bool polynomial::IsQuadratfrei() const -{ - polynomial Df(Derivative()); - if (Df.Degree() <= 0) { - return true; - } - if (GcdWith(Df).Degree() <= 1) { - return true; - } - else { - return false; - } -} - -template bool polynomial::CannotHaveRootsIn(const gInterval &I) const -{ - // get rid of easy cases - if (Degree() == -1) { - return false; - } - else if (Degree() == 0) { - return true; - } - else if (EvaluationAt(I.LowerBound()) == (T)0) { - return false; - } - - // generate list of derivatives - Gambit::List> derivs; - derivs.push_back(Derivative()); - int i; - for (i = 2; i <= Degree(); i++) { - derivs.push_back(derivs[i - 1].Derivative()); - } - - T val = EvaluationAt(I.LowerBound()); - if (val < (T)0) { - val = -val; - } - - // find max |c_0/Degree()*c_i|^(1/i) - int max_index = 0; - T base_of_max_index = (T)0; - T relevant_factorial = (T)1; - for (i = 1; i <= Degree(); i++) { - relevant_factorial *= (T)i; - T ith_coeff = derivs[i].EvaluationAt(I.LowerBound()) / relevant_factorial; - if (ith_coeff < (T)0) { - ith_coeff = -ith_coeff; - } - - if (ith_coeff != (T)0) { - T base = val / ((T)Degree() * ith_coeff); - - if (base_of_max_index == (T)0) { - max_index = i; - base_of_max_index = base; - } - else if (pow((T)base, max_index) < pow((T)base_of_max_index, i)) { - max_index = i; - base_of_max_index = base; - } - } - } - // assert(base_of_max_index != (T)0); - - if ((T)pow((T)I.Length(), max_index) < (T)base_of_max_index) { - return true; - } - else { - return false; - } -} - -template -Gambit::List> polynomial::RootSubintervals(const gInterval &I) const -{ // assert ( Degree() >= 0 && IsQuadratfrei() ); - - polynomial Df = Derivative(); - - Gambit::List> answer; - - Gambit::List> to_be_processed; - to_be_processed.push_back(I); - while (!to_be_processed.empty()) { - gInterval in_process = to_be_processed.back(); - to_be_processed.pop_back(); - - if (EvaluationAt(in_process.LowerBound()) == (T)0) { - if (Df.CannotHaveRootsIn(in_process)) { - answer.push_back(in_process); - } - else { - to_be_processed.push_back(in_process.RightHalf()); - to_be_processed.push_back(in_process.LeftHalf()); - } - } - else if (EvaluationAt(in_process.UpperBound()) == (T)0) { - if (Df.CannotHaveRootsIn(in_process)) { - if (in_process.UpperBound() == I.UpperBound()) { - answer.push_back(in_process); - } - } - else { - to_be_processed.push_back(in_process.RightHalf()); - to_be_processed.push_back(in_process.LeftHalf()); - } - } - - else if (!CannotHaveRootsIn(in_process)) { - if (Df.CannotHaveRootsIn(in_process)) { - if ((EvaluationAt(in_process.LowerBound()) < (T)0 && - EvaluationAt(in_process.UpperBound()) > (T)0) || - (EvaluationAt(in_process.LowerBound()) > (T)0 && - EvaluationAt(in_process.UpperBound()) < (T)0)) { - answer.push_back(in_process); - } - } - else { - to_be_processed.push_back(in_process.RightHalf()); - to_be_processed.push_back(in_process.LeftHalf()); - } - } - } - - return answer; -} - -template -gInterval polynomial::NeighborhoodOfRoot(const gInterval &I, T &error) const -{ - Gambit::List> intrvls; - intrvls.push_back(I); - while (intrvls.back().Length() >= error) { - if (EvaluationAt(intrvls.back().LowerBound()) == (T)0) { - intrvls.push_back(gInterval(intrvls.back().LowerBound(), intrvls.back().LowerBound())); - } - - else if (EvaluationAt(intrvls.back().UpperBound()) == (T)0) { - intrvls.push_back(gInterval(intrvls.back().UpperBound(), intrvls.back().UpperBound())); - } - - else if ((EvaluationAt(intrvls.back().LowerBound()) >= (T)0 && - EvaluationAt(intrvls.back().Midpoint()) <= (T)0) || - (EvaluationAt(intrvls.back().LowerBound()) <= (T)0 && - EvaluationAt(intrvls.back().Midpoint()) >= (T)0)) { - intrvls.push_back(intrvls.back().LeftHalf()); - } - - else { - intrvls.push_back(intrvls.back().RightHalf()); - } - } - - return intrvls.back(); - - // It is, perhaps, possible to speed this up, at least for double's - // by introducing Newton's method. -} - -template -Gambit::List> polynomial::PreciseRootIntervals(const gInterval &I, - T &error) const -{ - Gambit::List> coarse = RootSubintervals(I); - Gambit::List> fine; - - for (int i = 1; i <= coarse.size(); i++) { - fine.push_back(NeighborhoodOfRoot(coarse[i], error)); - } - - return fine; -} - -template -Gambit::List polynomial::PreciseRoots(const gInterval &I, T &error) const -{ - Gambit::List roots; - - polynomial p(*this), factor(*this); - - while (p.Degree() > 0) { - - int depth = 1; - polynomial probe(p.Derivative()); - polynomial current_gcd(p.GcdWith(probe)); - while (current_gcd.Degree() > 0) { - depth++; - factor = current_gcd; - probe = probe.Derivative(); - current_gcd = current_gcd.GcdWith(probe); - } - - for (int i = 1; i <= depth; i++) { - p = p / factor; - } - Gambit::List> fine = factor.PreciseRootIntervals(I, error); - for (int j = 1; j <= fine.size(); j++) { - T approx = fine[j].LowerBound(); - for (int h = 1; h <= 2; h++) { - approx -= factor.EvaluationAt(approx) / - factor.Derivative().EvaluationAt(approx); // Newton's Method - } - roots.push_back(approx); - } - factor = p; - } - - return roots; -} - -//-------------------------------------------------------------------------- -// class: complexpoly -//-------------------------------------------------------------------------- - -//-------------------------------------------------------------------------- -// constructors and a destructor -//-------------------------------------------------------------------------- - -complexpoly::complexpoly(const complexpoly &x) - - = default; - -complexpoly::complexpoly(const gComplex &coeff, const int °) -{ - if (coeff != (gComplex)0) { - for (int i = 0; i < deg; i++) { - coeflist.push_back((gComplex)0); - } - coeflist.push_back(coeff); - } -} - -complexpoly::complexpoly(const Gambit::List &coefficientlist) : coeflist(coefficientlist) -{ -} - -complexpoly::complexpoly(const int deg) : coeflist() -{ - if (deg >= 0) { - // gout << "Error is complexpoly int constructor.\n"; - exit(1); - } -} - -complexpoly::~complexpoly() = default; - -//-------------------------------------------------------------------------- -// operators -//-------------------------------------------------------------------------- - -complexpoly &complexpoly::operator=(const complexpoly &y) -{ - if (this != &y) { - coeflist = y.coeflist; - } - - return *this; -} - -bool complexpoly::operator==(const complexpoly &y) const -{ - if (Degree() != y.Degree()) { - return false; - } - else { - for (int i = 0; i <= Degree(); i++) { - if (coeflist[i + 1] != y.coeflist[i + 1]) { - return false; - } - } - } - return true; -} - -bool complexpoly::operator!=(const complexpoly &y) const { return !(*this == y); } - -const gComplex &complexpoly::operator[](const int index) const { return coeflist[index + 1]; } - -complexpoly complexpoly::operator+(const complexpoly &y) const -{ - if (Degree() < 0) { - return {y}; - } - else if (y.Degree() < 0) { - return {*this}; - } - - int max_degree; - - if (Degree() > y.Degree()) { - max_degree = Degree(); - } - else { - max_degree = y.Degree(); - } - - complexpoly sum; - for (int i = 0; i <= max_degree; i++) { - sum.coeflist.push_back((gComplex)0); - if (i <= Degree()) { - sum.coeflist[i + 1] += coeflist[i + 1]; - } - if (i <= y.Degree()) { - sum.coeflist[i + 1] += y.coeflist[i + 1]; - } - } - - while (!sum.coeflist.empty() && sum.coeflist.back() == (gComplex)0) { - sum.coeflist.pop_back(); - } - - return sum; -} - -complexpoly complexpoly::operator-(const complexpoly &y) const -{ - return complexpoly(*this + (-y)); -} - -complexpoly complexpoly::operator*(const complexpoly &y) const -{ - if (Degree() == -1) { - return {*this}; - } - else if (y.Degree() == -1) { - return {y}; - } - - int tot_degree = Degree() + y.Degree(); - - complexpoly product; - for (int t = 0; t <= tot_degree; t++) { - product.coeflist.push_back((gComplex)0); - } - for (int i = 0; i <= Degree(); i++) { - for (int j = 0; j <= y.Degree(); j++) { - product.coeflist[i + j + 1] += (*this)[i] * y[j]; - } - } - - return product; -} - -complexpoly complexpoly::operator/(const complexpoly &q) const -{ - // assert(q.Degree() >= 0); - - complexpoly ans; - complexpoly r = *this; - while (r.Degree() >= q.Degree()) { - complexpoly x(r.LeadingCoefficient() / q.LeadingCoefficient(), r.Degree() - q.Degree()); - ans += x; - r -= q * x; - } - return {ans}; -} - -complexpoly &complexpoly::operator+=(const complexpoly &y) { return ((*this) = (*this) + y); } - -complexpoly &complexpoly::operator-=(const complexpoly &y) { return ((*this) = (*this) - y); } - -complexpoly &complexpoly::operator*=(const complexpoly &y) { return ((*this) = (*this) * y); } - -complexpoly &complexpoly::operator/=(const complexpoly &y) { return ((*this) = (*this) / y); } - -complexpoly complexpoly::operator%(const complexpoly &q) const -{ - // assert (q.Degree() != -1); - - complexpoly ans; - complexpoly r = *this; - - while (r.Degree() >= q.Degree()) { - complexpoly x(r.LeadingCoefficient() / q.LeadingCoefficient(), r.Degree() - q.Degree()); - ans += x; - r -= q * x; - } - return {r}; -} - -complexpoly complexpoly::operator-() const -{ - complexpoly negation; - for (int i = 0; i <= Degree(); i++) { - negation.coeflist.push_back(-coeflist[i + 1]); - } - - return negation; -} - -complexpoly complexpoly::Derivative() const -{ - if (Degree() <= 0) { - return complexpoly(-1); - } - - complexpoly derivative; - for (int i = 1; i <= Degree(); i++) { - derivative.coeflist.push_back((gComplex)i * coeflist[i + 1]); - } - - return derivative; -} - -//-------------------------------------------------------------------------- -// manipulation -//-------------------------------------------------------------------------- - -void complexpoly::ToMonic() -{ - // assert (!IsZero()); - - gComplex lc = LeadingCoefficient(); - for (int i = 1; i <= coeflist.size(); i++) { - coeflist[i] /= lc; - } -} - -//-------------------------------------------------------------------------- -// information -//-------------------------------------------------------------------------- - -bool complexpoly::IsZero() const { return coeflist.empty(); } - -gComplex complexpoly::EvaluationAt(const gComplex &arg) const -{ - gComplex answer; - - for (int i = 0; i <= Degree(); i++) { - auto monom_val = static_cast(1); - for (int j = 1; j <= i; j++) { - monom_val *= arg; - } - answer += coeflist[i + 1] * monom_val; - } - - return answer; -} - -int complexpoly::Degree() const { return coeflist.size() - 1; } - -gComplex complexpoly::LeadingCoefficient() const -{ - if (Degree() < 0) { - return (gComplex)0; - } - else { - return coeflist[Degree() + 1]; - } -} - -complexpoly complexpoly::GcdWith(const complexpoly &that) const -{ - // assert( !this->IsZero() && !that.IsZero() ); - - complexpoly numerator(*this); - numerator.ToMonic(); - complexpoly denominator(that); - denominator.ToMonic(); - complexpoly remainder(numerator % denominator); - - while (!remainder.IsZero()) { - remainder.ToMonic(); - numerator = denominator; - denominator = remainder; - remainder = numerator % denominator; - } - - return denominator; -} - -bool complexpoly::IsQuadratfrei() const -{ - complexpoly Df(Derivative()); - if (Df.Degree() <= 0) { - return true; - } - if (GcdWith(Df).Degree() <= 1) { - return true; - } - else { - return false; - } -} - -Gambit::List complexpoly::Roots() const -{ - // assert (!IsZero()); - - Gambit::List answer; - - if (Degree() == 0) { - return answer; - } - - complexpoly deriv(Derivative()); - - gComplex guess(1.3, 0.314159); - - while (fabs(EvaluationAt(guess)) > 0.00001) { - gComplex diff = EvaluationAt(guess) / deriv.EvaluationAt(guess); - int count = 0; - bool done = false; - while (!done) { - if (count < 10 && fabs(EvaluationAt(guess - diff)) >= fabs(EvaluationAt(guess))) { - diff /= gComplex(4.0, 0.0); - count++; - } - else { - done = true; - } - } - if (count == 10) { - // gout << "Failure in complexpoly::Roots().\n"; - exit(1); - } - - guess -= diff; - } - - answer.push_back(guess); - - Gambit::List lin_form_coeffs; - lin_form_coeffs.push_back(guess); - lin_form_coeffs.push_back(gComplex(-1.0, 0.0)); - complexpoly linear_form(lin_form_coeffs); - complexpoly quotient = *this / linear_form; - for (const auto &root : quotient.Roots()) { - answer.push_back(root); - } - - for (int i = 1; i <= answer.size(); i++) { // "Polish" each root twice - answer[i] -= EvaluationAt(answer[i]) / deriv.EvaluationAt(answer[i]); - answer[i] -= EvaluationAt(answer[i]) / deriv.EvaluationAt(answer[i]); - } - - return answer; -} diff --git a/src/solvers/enumpoly/prepoly.cc b/src/solvers/enumpoly/prepoly.cc deleted file mode 100644 index 69b25e580..000000000 --- a/src/solvers/enumpoly/prepoly.cc +++ /dev/null @@ -1,497 +0,0 @@ -// -// This file is part of Gambit -// Copyright (c) 1994-2024, The Gambit Project (http://www.gambit-project.org) -// -// FILE: src/tools/enumpoly/prepoly.cc -// Implementation of supporting classes for polynomials -// -// 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 "gambit.h" - -#include "prepoly.h" - -//----------------------------------------------------------- -// gSpace -//----------------------------------------------------------- - -//------------------------- -// Constructors/Destructors -//------------------------- - -gSpace::gSpace(int nvars) : Variables() -{ - Variable *newvar; - for (int i = 1; i <= nvars; i++) { - newvar = new Variable; - newvar->Name = 'n'; - newvar->Name += Gambit::lexical_cast(i); - newvar->number = i; - Variables.push_back(newvar); - } -} - -gSpace::gSpace(const gSpace &p) : Variables() -{ - Variable *newvar; - for (int i = 1; i <= p.Variables.Length(); i++) { - newvar = new Variable; - newvar->Name = p.Variables[i]->Name; - newvar->number = i; - Variables.push_back(newvar); - } -} - -gSpace::~gSpace() -{ - for (int i = 1; i <= Variables.Length(); i++) { - delete Variables[i]; - } -} - -//----------------- -// Member Functions -//----------------- - -gSpace &gSpace::operator=(const gSpace &rhs) -{ - if (*this == rhs) { - return *this; - } - - Variables = rhs.Variables; - return *this; -} - -int gSpace::Dmnsn() const { return Variables.Length(); } - -Variable *gSpace::VariableWithNumber(int i) const { return Variables[i]; } - -Variable *gSpace::operator[](int i) const { return VariableWithNumber(i); } - -bool gSpace::operator==(const gSpace &rhs) const -{ - if (Variables.Length() == rhs.Variables.Length() && Variables == rhs.Variables) { - return true; - } - else { - return false; - } -} - -bool gSpace::operator!=(const gSpace &rhs) const { return !(*this == rhs); } - -//------------------------------------------------------ -// exp_vect -//------------------------------------------------------ - -//------------------------- -// Constructors/Destructors -//------------------------- - -exp_vect::exp_vect(const gSpace *p) : Space(p), components(p->Dmnsn()) -{ - for (int i = 1; i <= p->Dmnsn(); i++) { - components[i] = 0; - } -} - -exp_vect::exp_vect(const gSpace *p, const int &var, const int &exp) - : Space(p), components(p->Dmnsn()) -{ - for (int i = 1; i <= Dmnsn(); i++) { - components[i] = 0; - } - components[var] = exp; -} - -exp_vect::exp_vect(const gSpace *p, int *exponents) : Space(p), components(p->Dmnsn()) -{ - for (int i = 1; i <= Dmnsn(); i++) { - components[i] = exponents[i - 1]; - } -} - -exp_vect::exp_vect(const gSpace *p, Gambit::Vector exponents) - : Space(p), components(p->Dmnsn()) -{ - for (int i = 1; i <= Dmnsn(); i++) { - components[i] = exponents[i]; - } -} - -exp_vect::exp_vect(const gSpace *p, Gambit::Array exponents) - : Space(p), components(p->Dmnsn()) -{ - for (int i = 1; i <= Dmnsn(); i++) { - components[i] = exponents[i]; - } -} - -//------------------------- -// Operators -//------------------------- - -exp_vect &exp_vect::operator=(const exp_vect &RHS) -{ - if (this == &RHS) { - return *this; - } - - Space = RHS.Space; - components = RHS.components; - return *this; -} - -int exp_vect::operator[](int index) const { return components[index]; } - -bool exp_vect::operator==(const exp_vect &RHS) const -{ - if (components == RHS.components) { - return true; - } - else { - return false; - } -} - -bool exp_vect::operator!=(const exp_vect &RHS) const { return !(*this == RHS); } - -bool exp_vect::operator<=(const exp_vect &RHS) const -{ - for (int i = 1; i <= Dmnsn(); i++) { - if (components[i] > RHS.components[i]) { - return false; - } - } - - return true; -} - -bool exp_vect::operator>=(const exp_vect &RHS) const -{ - for (int i = 1; i <= Dmnsn(); i++) { - if (components[i] < RHS.components[i]) { - return false; - } - } - - return true; -} - -bool exp_vect::operator<(const exp_vect &RHS) const { return !(*this >= RHS); } - -bool exp_vect::operator>(const exp_vect &RHS) const { return !(*this <= RHS); } - -exp_vect exp_vect::operator-() const -{ - exp_vect tmp(Space); - for (int i = 1; i <= Dmnsn(); i++) { - tmp.components[i] = -components[i]; - } - - return tmp; -} - -exp_vect exp_vect::operator+(const exp_vect &credit) const -{ - // assert (Space == credit.Space); - - exp_vect tmp(Space); - for (int i = 1; i <= Dmnsn(); i++) { - tmp.components[i] = components[i] + credit.components[i]; - } - - return tmp; -} - -exp_vect exp_vect::operator-(const exp_vect &debit) const -{ - // assert (Space == debit.Space); - - exp_vect tmp(Space); - for (int i = 1; i <= Dmnsn(); i++) { - tmp.components[i] = components[i] - debit.components[i]; - } - - return tmp; -} - -void exp_vect::operator+=(const exp_vect &credit) -{ - - // assert (Space == credit.Space); - - for (int i = 1; i <= Dmnsn(); i++) { - components[i] += credit.components[i]; - } -} - -void exp_vect::operator-=(const exp_vect &debit) -{ - // assert (Space == debit.Space); - - for (int i = 1; i <= Dmnsn(); i++) { - components[i] -= debit.components[i]; - } -} - -//---------------------------- -// Other Operations -//---------------------------- - -exp_vect exp_vect::LCM(const exp_vect &arg2) const -{ - // assert (Space == arg2.Space); - - exp_vect tmp(Space); - for (int i = 1; i <= Dmnsn(); i++) { - if (components[i] < arg2.components[i]) { - tmp.components[i] = arg2.components[i]; - } - else { - tmp.components[i] = components[i]; - } - } - - return tmp; -} - -exp_vect exp_vect::AfterZeroingOutExpOfVariable(int &varnumber) const -{ - exp_vect tmp(*this); - tmp.components[varnumber] = 0; - return tmp; -} - -exp_vect exp_vect::AfterDecrementingExpOfVariable(int &varnumber) const -{ - exp_vect tmp(*this); - tmp.components[varnumber]--; - return tmp; -} - -//-------------------------- -// Information -//-------------------------- - -int exp_vect::Dmnsn() const { return Space->Dmnsn(); } - -bool exp_vect::IsConstant() const -{ - for (int i = 1; i <= Dmnsn(); i++) { - if ((*this)[i] > 0) { - return false; - } - } - return true; -} - -bool exp_vect::IsMultiaffine() const -{ - for (int i = 1; i <= Dmnsn(); i++) { - if ((*this)[i] > 1) { - return false; - } - } - return true; -} - -int exp_vect::TotalDegree() const -{ - int exp_sum = 0; - for (int i = 1; i <= Dmnsn(); i++) { - exp_sum += (*this)[i]; - } - return exp_sum; -} - -//-------------------------- -// Manipulation -//-------------------------- - -void exp_vect::ToZero() -{ - for (int i = 1; i <= Dmnsn(); i++) { - components[i] = 0; - } -} - -//------------------------------------------------------ -// term_order -//------------------------------------------------------ - -//------------------- -// Possible Orderings -//------------------- - -bool lex(const exp_vect &LHS, const exp_vect &RHS) -{ - for (int i = 1; i <= LHS.Dmnsn(); i++) { - if (LHS[i] < RHS[i]) { - return true; - } - else if (LHS[i] > RHS[i]) { - return false; - } - } - return false; -} - -bool reverselex(const exp_vect &LHS, const exp_vect &RHS) -{ - for (int i = LHS.Dmnsn(); i >= 1; i--) { - if (LHS[i] < RHS[i]) { - return true; - } - else if (LHS[i] > RHS[i]) { - return false; - } - } - return false; -} - -bool deglex(const exp_vect &LHS, const exp_vect &RHS) -{ - if (LHS.TotalDegree() < RHS.TotalDegree()) { - return true; - } - else if (LHS.TotalDegree() > RHS.TotalDegree()) { - return false; - } - - for (int i = 1; i <= LHS.Dmnsn(); i++) { - if (LHS[i] < RHS[i]) { - return true; - } - else if (LHS[i] > RHS[i]) { - return false; - } - } - return false; -} - -bool reversedeglex(const exp_vect &LHS, const exp_vect &RHS) -{ - if (LHS.TotalDegree() < RHS.TotalDegree()) { - return true; - } - else if (LHS.TotalDegree() > RHS.TotalDegree()) { - return false; - } - - for (int i = LHS.Dmnsn(); i >= 1; i--) { - if (LHS[i] < RHS[i]) { - return true; - } - else if (LHS[i] > RHS[i]) { - return false; - } - } - return false; -} - -bool degrevlex(const exp_vect &LHS, const exp_vect &RHS) -{ - if (LHS.TotalDegree() < RHS.TotalDegree()) { - return true; - } - else if (LHS.TotalDegree() > RHS.TotalDegree()) { - return false; - } - - for (int i = LHS.Dmnsn(); i >= 1; i--) { - if (LHS[i] < RHS[i]) { - return false; - } - else if (LHS[i] > RHS[i]) { - return true; - } - } - return false; -} - -bool reversedegrevlex(const exp_vect &LHS, const exp_vect &RHS) -{ - if (LHS.TotalDegree() < RHS.TotalDegree()) { - return true; - } - else if (LHS.TotalDegree() > RHS.TotalDegree()) { - return false; - } - - for (int i = 1; i <= LHS.Dmnsn(); i++) { - if (LHS[i] < RHS[i]) { - return false; - } - else if (LHS[i] > RHS[i]) { - return true; - } - } - return false; -} - -//------------------------- -// Constructors/Destructors -//------------------------- - -term_order::term_order(const gSpace *p, ORD_PTR act_ord) : Space(p), actual_order(act_ord) {} - -//------------------------- -// Operators -//------------------------- - -term_order &term_order::operator=(const term_order &RHS) -{ - if (this == &RHS) { - return *this; - } - - Space = RHS.Space; - actual_order = RHS.actual_order; - return *this; -} - -bool term_order::operator==(const term_order &RHS) const -{ - return (Space == RHS.Space && actual_order == RHS.actual_order); -} - -bool term_order::operator!=(const term_order &RHS) const { return !(*this == RHS); } - -//------------------------- -// Comparisons -//------------------------- - -bool term_order::Less(const exp_vect &LHS, const exp_vect &RHS) const -{ - return (*actual_order)(LHS, RHS); -} - -bool term_order::LessOrEqual(const exp_vect &LHS, const exp_vect &RHS) const -{ - return ((*actual_order)(LHS, RHS) || LHS == RHS); -} - -bool term_order::Greater(const exp_vect &LHS, const exp_vect &RHS) const -{ - return !(LessOrEqual(LHS, RHS)); -} - -bool term_order::GreaterOrEqual(const exp_vect &LHS, const exp_vect &RHS) const -{ - return !(Less(LHS, RHS)); -} diff --git a/src/solvers/enumpoly/prepoly.h b/src/solvers/enumpoly/prepoly.h deleted file mode 100644 index 18be0074a..000000000 --- a/src/solvers/enumpoly/prepoly.h +++ /dev/null @@ -1,182 +0,0 @@ -// -// This file is part of Gambit -// Copyright (c) 1994-2024, The Gambit Project (http://www.gambit-project.org) -// -// FILE: src/tools/enumpoly/prepoly.h -// Declaration of supporting classes for polynomials -// -// 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. -// - -#ifndef PREPOLY_H -#define PREPOLY_H - -#include -#include "gambit.h" - -/* - The classes in this file are prior to the notion of a -multivariate polynomial. First of all, one needs a space which the -polynomials refer to. This is given by the notion of a gSpace. -A polynomial is a sum of monomials, where each monomial is a -coefficient multiplied by a power product, and each power product is -the vector of variables "raised" by an exponent vector. The exponent -vectors play a special role in many algorithms, especially those -associated with Groebner bases, and they have a class. Finally, the -notion of a term order on the set of exponent vector, which -encompasses potentially infinitely many distinct orders, seems best -implemented through a class that delivers a full menu of services -built on top of a pointer to a function for computing an order. -*/ - -// ************************* -// gSpace declaration -// ************************* - -struct Variable { - std::string Name; - int number; -}; - -class gSpace { -private: - Gambit::Array Variables; - -public: - explicit gSpace(int nvars = 0); - gSpace(const gSpace &); - ~gSpace(); - - // operators - gSpace &operator=(const gSpace &rhs); - - Variable *operator[](int) const; - bool operator==(const gSpace &rhs) const; - bool operator!=(const gSpace &rhs) const; - - // information - int Dmnsn() const; - Variable *VariableWithNumber(int) const; -}; - -// *********************** -// class exp_vect -// *********************** - -/* - Exponent vectors are vectors of integers. In specifying operators - we take the view that addition, subtraction, and order are defined, - but not inner or scalar products. -*/ - -class exp_vect { - -private: - const gSpace *Space; - Gambit::Vector components; - -public: - explicit exp_vect(const gSpace *); - exp_vect(const gSpace *, const int &, const int &); // x_i^j - exp_vect(const gSpace *, int *); - exp_vect(const gSpace *, Gambit::Vector); - exp_vect(const gSpace *, Gambit::Array); - exp_vect(const exp_vect &) = default; - ~exp_vect() = default; - - // Operators - exp_vect &operator=(const exp_vect &RHS); - - int operator[](int index) const; - - bool operator==(const exp_vect &RHS) const; - bool operator!=(const exp_vect &RHS) const; - bool operator<=(const exp_vect &RHS) const; - bool operator>=(const exp_vect &RHS) const; - bool operator<(const exp_vect &RHS) const; - bool operator>(const exp_vect &RHS) const; - - exp_vect operator-() const; - exp_vect operator+(const exp_vect &) const; - exp_vect operator-(const exp_vect &) const; - void operator+=(const exp_vect &); - void operator-=(const exp_vect &); - - // Other operations - exp_vect LCM(const exp_vect &) const; - exp_vect AfterZeroingOutExpOfVariable(int &) const; - exp_vect AfterDecrementingExpOfVariable(int &) const; - - // Information - int Dmnsn() const; - bool IsConstant() const; - bool IsMultiaffine() const; - int TotalDegree() const; - - // Manipulation - void ToZero(); -}; - -// *********************** -// class term_order -// *********************** - -/* - A term order is a total order of the set of exponent vectors -associated with a particular variable list, that has the properties: - - a) 1 < alpha for all alpha \ne 1; - b) if alpha < beta, then alpha + gamma < beta + gamma for all gamma >= 0. - - In our implementation we take the view that the order itself is a -variable of an object of the class, and implement this in terms of -pointers to functions. -*/ - -// THE FOLLOWING FUNCTIONS SHOULD BE VIEWED AS PRIVATE MEMBERS OF -// class term_order I WAS BAFFLED AS TO HOW TO HAVE A MEMBER THAT -// IS A POINTER-TO-OTHER-MEMBER-FUNCTION -typedef bool (*ORD_PTR)(const exp_vect &, const exp_vect &); -bool lex(const exp_vect &, const exp_vect &); -bool reverselex(const exp_vect &, const exp_vect &); -bool deglex(const exp_vect &, const exp_vect &); -bool reversedeglex(const exp_vect &, const exp_vect &); -bool degrevlex(const exp_vect &, const exp_vect &); -bool reversedegrevlex(const exp_vect &, const exp_vect &); - -class term_order { -private: - const gSpace *Space; - ORD_PTR actual_order; - -public: - term_order(const gSpace *, ORD_PTR); - term_order(const term_order &) = default; - ~term_order() = default; - - // Operators - term_order &operator=(const term_order &RHS); - - bool operator==(const term_order &RHS) const; - bool operator!=(const term_order &RHS) const; - - // Comparisons invoking the underlying order - bool Less(const exp_vect &, const exp_vect &) const; - bool LessOrEqual(const exp_vect &, const exp_vect &) const; - bool Greater(const exp_vect &, const exp_vect &) const; - bool GreaterOrEqual(const exp_vect &, const exp_vect &) const; -}; - -#endif // PREPOLY_H diff --git a/src/solvers/enumpoly/quiksolv.h b/src/solvers/enumpoly/quiksolv.h index 292f00340..a6ec4515a 100644 --- a/src/solvers/enumpoly/quiksolv.h +++ b/src/solvers/enumpoly/quiksolv.h @@ -123,8 +123,7 @@ template class QuikSolv { bool operator!=(const QuikSolv &) const; // Information - inline const gSpace *AmbientSpace() const { return System.AmbientSpace(); } - inline const term_order *TermOrder() const { return System.TermOrder(); } + inline const VariableSpace *AmbientSpace() const { return System.AmbientSpace(); } inline int Dmnsn() const { return System.Dmnsn(); } inline const gPolyList &UnderlyingEquations() const { return System; } inline bool WasSolved() const { return HasBeenSolved; } diff --git a/src/solvers/enumpoly/quiksolv.imp b/src/solvers/enumpoly/quiksolv.imp index 8154757d2..e13beb0d6 100644 --- a/src/solvers/enumpoly/quiksolv.imp +++ b/src/solvers/enumpoly/quiksolv.imp @@ -34,7 +34,7 @@ using namespace Gambit; template QuikSolv::QuikSolv(const gPolyList &given) - : System(given), gDoubleSystem(given.AmbientSpace(), given.TermOrder(), given.NormalizedList()), + : System(given), gDoubleSystem(given.AmbientSpace(), given.NormalizedList()), NoEquations(std::min(System.Dmnsn(), System.Length())), NoInequalities(std::max(System.Length() - System.Dmnsn(), 0)), TreesOfPartials(gDoubleSystem), HasBeenSolved(false), Roots(), isMultiaffine(System.IsMultiaffine()),