gpoly.imp

Gambit 是一个游戏库理论软件

	else activar = j;
    }
  }
  return activar;
}

template <class T> 
polynomial<T> gPoly<T>::UnivariateEquivalent(int activar) const
{
  assert(UniqueActiveVariable() >= 0);

  gList<T> coefs;

  if (!IsZero()) {
    for (int h = 0; h <= DegreeOfVar(activar); h++) coefs += (T)0;

    for (int i = 1; i <= Terms.Length(); i++)
      coefs[Terms[i].ExpV()[activar] + 1] = Terms[i].Coef();
  }

  return  polynomial<T>(coefs);
} 

//-------------------------------------------------------------
//           Private Versions of Arithmetic Operators
//-------------------------------------------------------------

template <class T> 
gList< gMono<T> > gPoly<T>::Adder(const gList<gMono<T> >& One, 
				   const gList<gMono<T> >& Two) const
{
  if (One.Length() == 0) return Two; if (Two.Length() == 0) return One;

  gList<gMono<T> > answer;

  int i = 1; int j = 1;
  while (i <= One.Length() || j <= Two.Length()) {
    if      (i > One.Length()) { answer += Two[j]; j++; }
    else if (j > Two.Length()) { answer += One[i]; i++; }
    else {
      if      ( Order->Less(   One[i].ExpV(),Two[j].ExpV()) ) 
	{ answer += One[i]; i++; }
      else if ( Order->Greater(One[i].ExpV(),Two[j].ExpV()) ) 
	{ answer += Two[j]; j++; }
      else {
	if (One[i].Coef() + Two[j].Coef() != (T)0) 
	  answer += One[i] + Two[j]; 
	i++; j++; 
      }
    }
  }

  return answer;
}

template <class T> 
gList< gMono<T> > gPoly<T>::Mult(const gList<gMono<T> >& One, 
				  const gList<gMono<T> >& Two) const
{
  gList<gMono<T> > answer;

  if (One.Length() == 0 || Two.Length() == 0) return answer;

  int i;
  for (i = 1; i <= One.Length(); i++)
    for (int j = 1; j <= Two.Length(); j++) 
      {
	gMono<T> next = One[i] * Two[j];

	if (answer.Length() == 0) answer += next;

	else
	  {
	    int bot = 1; int top = answer.Length();
	    if      ( Order->Less   (answer[top].ExpV(),next.ExpV()) ) 
	      answer += next;
	    else if ( Order->Greater(answer[bot].ExpV(),next.ExpV()) )
	      answer.Insert(next,1);
	    else
	      {
		if      ( answer[bot].ExpV() == next.ExpV() ) top = bot;
		else if ( answer[top].ExpV() == next.ExpV() ) bot = top;
		
		while (bot < top - 1)
		  {
		    int test = bot + (top - bot)/2;
		    if ( answer[test].ExpV() == next.ExpV() ) bot = top = test;
		    else
		      if      (Order->Less   (answer[test].ExpV(),next.ExpV()))
			bot = test;
		      else // (Order->Greater(answer[test].ExpV(),next.ExpV()))
			top = test;
		  }

		if (bot == top) answer[bot] += next;
		else            answer.Insert(next,top);
	      }
	  }
      }
  return answer;
}

template<class T> 
gPoly<T> gPoly<T>::DivideByPolynomial(const gPoly<T> &den) const
{
  gPoly<T> zero(Space,(T)0,Order);

  assert(den != zero);
  assert(*this == zero || den.Degree() <= Degree());

  // assumes exact divisibility!

  gPoly<T> result = zero;

  if (*this == zero) return result;
  else if (den.Degree() == 0)
    { result = *this/den.NumLeadCoeff(); return result; }
  else
    {
      int last = Dmnsn();
      while (den.DegreeOfVar(last) == 0) last--;

      gPoly<T> remainder = *this;

      while (remainder != zero)
	{
	  gPoly<T> quot = remainder.LeadingCoefficient(last) /
	                     den.LeadingCoefficient(last);
	  gPoly<T> power_of_last(Space,last,remainder.DegreeOfVar(last) - 
			                        den.DegreeOfVar(last),Order);
	  result += quot * power_of_last;
	  remainder -= quot * power_of_last * den;
	}
    }
  return result;
} 

template <class T> 
gPoly<T> gPoly<T>::EvaluateOneVar( int varnumber, T val) const
{
  gPoly<T> answer(Space,(T)0,Order);

  for (int i = 1; i <= Terms.Length(); i++) 
    answer += gPoly<T>(Space,
	       Terms[i].ExpV().AfterZeroingOutExpOfVariable(varnumber),
	       ((T) Terms[i].Coef())* ((T) pow(val,(long int)Terms[i].ExpV()[varnumber])), 
		       Order);
  return answer;
}

template <class T>
exp_vect gPoly<T>::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.Length(); i++)
    if ( order.Less(answer,Terms[i].ExpV()) && answer < Terms[i].ExpV() )
      answer = Terms[i].ExpV();
  return answer;
}

template <class T> 
gPoly<T> gPoly<T>::PartialDerivative(int varnumber) const
{
  gPoly<T> newPoly(*this); 

  for (int i = 1; i <= newPoly.Terms.Length(); i++)
    newPoly.Terms[i] = gMono<T>(newPoly.Terms[i].Coef()
                           * (T)newPoly.Terms[i].ExpV()[varnumber],
	  newPoly.Terms[i].ExpV().AfterDecrementingExpOfVariable(varnumber));

  int j = 1;
  while (j <= newPoly.Terms.Length()) {
    if (newPoly.Terms[j].Coef() == (T)0) newPoly.Terms.Remove(j);
    else j++;
  }

  return (newPoly);
} 

template <class T> 
gPoly<T> gPoly<T>::LeadingCoefficient(int varnumber) const
{
  gPoly<T> newPoly(*this); 

  int degree = DegreeOfVar(varnumber);

  newPoly.Terms = gList<gMono<T> >();
  for (int j = 1; j <= Terms.Length(); j++) {
    if (Terms[j].ExpV()[varnumber] == degree)
      newPoly.Terms += gMono<T>(Terms[j].Coef(),
	      Terms[j].ExpV().AfterZeroingOutExpOfVariable(varnumber));
  }

  return (newPoly);
}

//--------------------
// Term Order Concepts
//--------------------

template <class T> 
exp_vect gPoly<T>::LeadingPowerProduct(const term_order & order) const
{
  assert (Terms.Length() > 0);

  if (*Order == order) // worth a try ...
    return Terms[Terms.Length()].ExpV();
  else {
    int max = 1;
    for (int j = 2; j <= Terms.Length(); j++) {
      if ( order.Less(Terms[max].ExpV(),Terms[j].ExpV()) )
	max = j;
    }
    return Terms[max].ExpV();
  }
}

template <class T> 
T gPoly<T>::LeadingCoefficient(const term_order & order) const
{
  if (*Order == order) // worth a try ...
    return Terms[Terms.Length()].Coef();
  else {
    int max = 1;
    for (int j = 2; j <= Terms.Length(); j++)
      if ( order.Less(Terms[max].ExpV(),Terms[j].ExpV()) )
	max = j;
    return Terms[max].Coef();
  }
}

template <class T>
gPoly<T> gPoly<T>::LeadingTerm(const term_order & order) const
{
  if (*Order == order) // worth a try ...
    return gPoly<T>(Space,Terms[Terms.Length()],Order);
  else {
    int max = 1;
    for (int j = 2; j <= Terms.Length(); j++)
      if ( order.Less(Terms[max].ExpV(),Terms[j].ExpV()) )
	max = j;
    return gPoly<T>(Space,Terms[max],Order);
  }
}

template <class T>
void gPoly<T>::ToMonic(const term_order & order) 
{
  *this = *this/LeadingCoefficient(order);
}

template <class T>
void gPoly<T>::ReduceByDivisionAtExpV(const term_order & order, 
				      const gPoly<T> & divisor, 
				      const exp_vect & expv)
{
  assert(expv >= divisor.LeadingPowerProduct(order));
  assert(divisor.LeadingCoefficient(order) != (T)0);

  gPoly<T> factor(Space, 
		    expv - divisor.LeadingPowerProduct(order), 
		    (T)1,
		    Order);

  *this -= (GetCoef(expv) / divisor.LeadingCoefficient(order)) * 
           factor * divisor;
}

template <class T>
void gPoly<T>::ReduceByRepeatedDivision(const term_order & order, 
					const gPoly<T> & 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 <class T>
gPoly<T> gPoly<T>::S_Polynomial(const term_order& order, 
                                const gPoly<T>& arg2) const
{
  exp_vect exp_lcm = 
        (LeadingPowerProduct(order)).LCM(arg2.LeadingPowerProduct(order));
  gPoly<T> lcm = gPoly<T>(Space,exp_lcm,(T)1,Order);
  gPoly<T> L1 = lcm/LeadingTerm(order);
  gPoly<T> L2 = lcm/arg2.LeadingTerm(order);

  return L1*(*this) - L2*arg2;
}

template<class T> gPoly<T> 
gPoly<T>::TranslateOfMono(const gMono<T>& m, 
				      const gVector<T>& new_origin) const
{
  gPoly<T> answer(GetSpace(), (T)1, GetOrder());

  for (int i = 1; i <= Dmnsn(); i++) {
    if (m.ExpV()[i] > 0) {
      gPoly<T> lt(GetSpace(), i, 1, GetOrder());
      lt += gPoly<T>(GetSpace(), new_origin[i], GetOrder());
      for (int j = 1; j <= m.ExpV()[i]; j++)
	answer *= lt;
    }
  }

  answer *= m.Coef();

  return answer;
}

template<class T> gPoly<T> 
gPoly<T>::TranslateOfPoly(const gVector<T>& new_origin) const
{
  gPoly<T> answer(GetSpace(), GetOrder());
  for (int i = 1; i <= this->MonomialList().Length(); i++) 
    answer += TranslateOfMono(this->MonomialList()[i],new_origin);
  return answer;
}


template<class T> gPoly<T>  
gPoly<T>::MonoInNewCoordinates(const gMono<T>& m, 
					    const gSquareMatrix<T>& M) const
{
  assert(M.NumRows() == Dmnsn());

  gPoly<T> answer(GetSpace(), (T)1, GetOrder());

  for (int i = 1; i <= Dmnsn(); i++) {
    if (m.ExpV()[i] > 0) {
      gPoly<T> linearform(GetSpace(), (T)0, GetOrder());
      for (int j = 1; j <= Dmnsn(); j++) {
	exp_vect exps(GetSpace(), j, 1);
	linearform += gPoly<T>(GetSpace(), exps, M(i,j), GetOrder());
      }
      for (int k = 1; k <= m.ExpV()[i]; k++) answer *= linearform;
    }
  }

  answer *= m.Coef();

  return answer;
}

template<class T> 
gPoly<T> gPoly<T>::PolyInNewCoordinates(const gSquareMatrix<T>& M) const
{
  gPoly<T> answer(GetSpace(), GetOrder());
  for (int i = 1; i <= MonomialList().Length(); i++) 
    answer += MonoInNewCoordinates(MonomialList()[i],M);
  return answer;
}

template<class T>
T gPoly<T>::MaximalValueOfNonlinearPart(const T& radius)  const
{
  T maxcon = (T)0;
  for (int i = 1; i <= MonomialList().Length(); i++) 
    if (MonomialList()[i].TotalDegree() > 1) 
      maxcon += MonomialList()[i].Coef() * 
	            pow(radius,(long)MonomialList()[i].TotalDegree());

  return maxcon;
}

//---------------------------
//    Global Operators
//---------------------------

template<class T> gPoly<T> operator*(const T val, const gPoly<T> &poly)
{
  gPoly<T> prod(poly);
  prod *= val;
  return prod;
}

template<class T> gPoly<T> operator*(const gPoly<T> &poly, const T val)
{
  return val * poly;
}

template<class T> gPoly<T> operator+(const T val, const gPoly<T> &poly)
{
  gPoly<T> prod(poly);
  prod += val;
  return prod;
}

template<class T> gPoly<T> operator+(const gPoly<T> &poly, const T val)
{
  return val + poly;
}

template <class T> void gPoly<T>::Output(gText &t) const
{
  gText s;
  if (Terms.Length() == 0)  {
    s += "0";
  }
  else {
    for (int j = 1; j <= Terms.Length(); 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 += ToText(-Terms[j].Coef());
	}
	else {                        
	  s += ToText( 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 += ToText(exp);
	    }
	  }
	}
	
	if (j < Terms.Length()) s += " ";
      }
  }
  if (s == "") s = " 0";
  
  t += s;
}

template <class T> gText &operator<<(gText &p_text, const gPoly<T> &p_poly)
{
  p_poly.Output(p_text);
  return p_text;
}

//----------------------------------
//           Conversion
//----------------------------------

template<class T>  
gPoly<gDouble>  TogDouble(const gPoly<T>& given) 
{
  gPoly<gDouble>   answer(given.GetSpace(),given.GetOrder());
  gList<gMono<T> > list = given.MonomialList();
  for (int i =1; i <= list.Length(); i++) {
    gDouble nextcoef = (gDouble)list[i].Coef();
    gPoly<gDouble> next(given.GetSpace(),
			list[i].ExpV(),
			nextcoef,
			given.GetOrder());
    answer += next;
  }
  
  return answer;
}

template<class T>  
gPoly<gDouble>  NormalizationOfPoly(const gPoly<T>& given) 
{
  gList<gMono<T> > list = given.MonomialList();
  gDouble maxcoeff = (gDouble)0;
  for (int i =1; i <= list.Length(); i++) {
    maxcoeff = gmax(maxcoeff,abs((gDouble)list[i].Coef()));
  }

  if (maxcoeff < 0.000000001) return TogDouble(given);

  gPoly<gDouble>   answer(given.GetSpace(),given.GetOrder());
  for (int i =1; i <= list.Length(); i++) {
    gDouble nextcoef = (gDouble)list[i].Coef()/maxcoeff;
    gPoly<gDouble> next(given.GetSpace(),
			list[i].ExpV(),
			nextcoef,
			given.GetOrder());
    answer += next;
  }
  
  return answer;
}

template <class T> gOutput &operator<<(gOutput &f, const gPoly<T> &p)
{
  gText s;
  s << p;
  f << s;
  return f;
}