gpoly.imp

Gambit 是一个游戏库理论软件

//
// $Source: /home/gambit/CVS/gambit/sources/poly/gpoly.imp,v $
// $Date: 2002/08/27 17:29:46 $
// $Revision: 1.3 $
//
// DESCRIPTION:
// Implementation of multivariate polynomial type
//
// This file is part of Gambit
// Copyright (c) 2002, The Gambit Project
//
// 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 "gpoly.h"
#include "math/gmath.h"

//---------------------------------------------------------------
//                      gPoly
//---------------------------------------------------------------

//---------------------------
// Constructors / Destructor
//---------------------------


template <class T> gPoly<T>::gPoly(const gSpace* p, const term_order* o)
: Space(p), Order(o), Terms()
{
}

template <class T> gPoly<T>::gPoly(const gSpace*     p, 
				     const T&          constant,
				     const term_order* o)
: Space(p), Order(o), Terms()
{
  if (constant != (T)0)   Terms += gMono<T>(p,constant);
}

template <class T> gPoly<T>::gPoly(const gSpace *p, 
				       const gText &s,
				       const term_order* o)
: Space(p), Order(o), Terms()
{
  *this=s;  // Operator = needs to be fixed
}

template <class T> gPoly<T>::gPoly(const gSpace *p, int var_no, int exp,
				       const term_order* o)
: Space(p), Order(o), Terms()
{
  Terms += gMono<T>((T)1,exp_vect(p,var_no,exp));
}

template <class T> gPoly<T>::gPoly(const gSpace* p, 
				       exp_vect exps, 
				       T coeff,
				       const term_order* o)
: Space(p), Order(o), Terms()
{
  Terms += gMono<T>(coeff,exps);
}

template <class T> gPoly<T>::gPoly(const gSpace* p, 
				   const gMono<T>& mono, 
				   const term_order* o)
: Space(p), Order(o), Terms()
{
  Terms += mono;
}

template <class T> gPoly<T>::gPoly(const gPoly<T> &p)
: Space(p.Space), Order(p.Order), Terms(p.Terms)
{
  *this=p; 
}

template <class T> gPoly<T>::~gPoly()
{
}


//----------------------------------
//        Operators
//----------------------------------

template<class T> gPoly<T> &gPoly<T>::operator=(const gPoly<T> &p)
{
  assert (Space == p.Space && Order == p.Order);

  Terms = p.Terms; return (*this);
}

template <class T> gPoly<T>& gPoly<T>::operator=(const gText &Hold)
{
  gList<gMono<T> > nullTerms;
  Terms = nullTerms;           // get rid of old Terms

  charnum = 0;
  int contflag = 1;
  T nega = 1;
  gArray<int> 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;
    }
  }

  gList<gMono<T> > newTerms;
  for (int j = 1; j <= Terms.Length(); j++) {
    int low = 0; int high = newTerms.Length() + 1;
    while (low +1 < high) {
      int test = low + (high - low)/2;
      if (1 <= test && test <= newTerms.Length()) 
	assert (Terms[j].ExpV() != newTerms[test].ExpV());
      if ( Order->Less(Terms[j].ExpV(),newTerms[test].ExpV()) )
	high = test;
      else
	low = test;
    }
    newTerms.Insert(Terms[j],high);
  }
  Terms = newTerms;

  return (*this);
}


template<class T> gPoly<T> gPoly<T>::operator-() const
{
  gPoly<T> neg(*this);
  for (int j = 1; j <= Terms.Length(); j++)  neg.Terms[j] = -Terms[j];
  return neg;
}

template<class T> gPoly<T> gPoly<T>::operator-(const gPoly<T> &p) const
{
  gPoly<T> dif(*this);
  dif -= p;
  return dif;
}

template<class T> void gPoly<T>::operator-=(const gPoly<T> &p)
{
  assert(Space == p.Space);

  gPoly<T> neg = p;
  for (int i = 1; i <= neg.Terms.Length(); i++)
    neg.Terms[i] = - neg.Terms[i];
  Terms = Adder(Terms,neg.Terms);
}

template<class T> gPoly<T> gPoly<T>::operator+(const gPoly<T> &p) const
{ 
  gPoly<T> sum(*this);
  sum += p;
  return sum;
}

template<class T> void gPoly<T>::operator+=(const gPoly<T> &p)
{
  assert(Space == p.Space);

  Terms = Adder(Terms,p.Terms); 
}

template<class T> void gPoly<T>::operator+=(const T& val)
{
  *this += gPoly<T>(Space,val,Order);
}

template<class T> gPoly<T> gPoly<T>::operator*(const gPoly<T> &p) const
{
  gPoly<T> prod(*this);
  prod *= p;
  return prod;
}

template<class T> gPoly<T> gPoly<T>::operator/(const T val) const
{
  assert (val != (T)0);
  T one = (T)1;
  return (one/val) * (*this);
}

template<class T> 
gPoly<T> gPoly<T>::operator/(const gPoly<T> & den) const
{
  return DivideByPolynomial(den);
}

template<class T> void gPoly<T>::operator*=(const gPoly<T> &p)
{ 
  assert(Space == p.Space);

  Terms = Mult(Terms,p.Terms);
}

template<class T> void gPoly<T>::operator*=(const T& val)
{ 
  for (int j = 1; j <= Terms.Length(); j++)
    Terms[j] *= val;
}


template<class T> bool gPoly<T>::operator==(const gPoly<T> &p) const
{
  assert(Space == p.Space && Order == p.Order);

  if (Terms.Length() != p.Terms.Length())           return false;
  if (Terms.Length() == 0 && p.Terms.Length() == 0) return true;

  return (Terms == p.Terms);
}

template<class T> bool gPoly<T>::operator!=(const gPoly<T> &p) const
{
  return !(*this == p);
}


//----------------------------------
//        Member Functions
//----------------------------------

template <class T> void gPoly<T>::String_Term(T nega)
{
  gArray<int> 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 += gMono<T>(val, exp_vect(Space,PowArray));
}

template <class T> int gPoly<T>::String_GetPow(void)
{

  gText Pow = "";
  while (charc == ' '){
    charnum++;
    charc = TheString[charnum];
  }

  while (charc >= '0' && charc <= '9'){
    Pow += charc;
    charnum++;
    charc = TheString[charnum];
  }
  return (atoi(Pow));
}

template <class T> void gPoly<T>::String_VarAndPow(gArray<int> &PowArray)
{
  gText 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 <class T> bool gPoly<T>::Check_String(const gText &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}; 
 gText 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 <class T> void gPoly<T>::GetChar(void)
{
  charc = TheString[charnum];
}

//----------------------------------
//           Information
//----------------------------------
  
template <class T> const gSpace* gPoly<T>::GetSpace(void) const
{
  return (gSpace *) Space;  
}
  
template <class T> const term_order* gPoly<T>::GetOrder(void) const
{
  return (term_order *) Order;  
}

template <class T> int gPoly<T>::Dmnsn() const
{
  return Space->Dmnsn();
}

template <class T> int gPoly<T>::DegreeOfVar(int var_no) const
{
  int max = 0;
  for (int j = 1; j <= Terms.Length(); j++) 
    if (max < Terms[j].ExpV()[var_no]) 
      max = Terms[j].ExpV()[var_no];
  return max;
}

template <class T> bool gPoly<T>::IsZero() const
{
  if (Terms.Length() == 0) return true;
  else                     return false;
}

template <class T> int gPoly<T>::Degree() const
{
  int max = 0;
  for (int j = 1; j <= Terms.Length(); j++)
    if (Terms[j].TotalDegree() > max)  
      max = Terms[j].TotalDegree();
  return max;
}

template <class T> T gPoly<T>::GetCoef(const gArray<int> &Powers) const
{
  return GetCoef(exp_vect(Space,Powers));
}

template <class T> T gPoly<T>::GetCoef(const exp_vect &Powers) const
{
  for (int j = 1; j <= Terms.Length(); j++)
    if (Terms[j].ExpV() == Powers)
      return Terms[j].Coef();
  return (T)0;
}

template <class T> T gPoly<T>::NumLeadCoeff() const
{
  assert (Degree() == 0 && Terms.Length() <= 1);

  if (Terms.Length() == 1) return Terms[1].Coef();
  else                     return (T)0;
}

template <class T> bool gPoly<T>::IsConstant() const
{
  for (int i = 1; i <= Terms.Length(); i++)
    if (!Terms[i].IsConstant())
      return false;
  return true;
}

template <class T> bool gPoly<T>::IsMultiaffine() const
{
  for (int i = 1; i <= Terms.Length(); i++)
    if (!Terms[i].IsMultiaffine())
      return false;
  return true;
}

template <class T> T gPoly<T>::Evaluate(const gArray<T> &values) const
{
  T answer=0;
  for (int j = 1; j <= Terms.Length(); j++)
    answer += Terms[j].Evaluate(values);
  return answer;
}


template <class T> gList<exp_vect> gPoly<T>::ExponentVectors() const
{
  gList<exp_vect> result;
  for (int j = 1; j <= Terms.Length(); j++)
    result += exp_vect(Terms[j].ExpV());
  return result;
}

template <class T> gList<gMono<T> > gPoly<T>::MonomialList() const
{
  return Terms;
}

template <class T> int gPoly<T>::No_Monomials() const
{
  //  gout << "Eliminate old code in No_monomials, if successful.\n";

  return Terms.Length();
}

template <class T> int gPoly<T>::UniqueActiveVariable() const
{
  gList<exp_vect> ExpVecs = ExponentVectors();
  int activar = 0;
  for (int i = 1; i <= ExpVecs.Length(); i++) {
    for (int j = 1; j <= Dmnsn(); j++) {
      if (ExpVecs[i][j] > 0)
	if (activar > 0 && activar != j) return -1;// multivariate!