poly.imp

Gambit 是一个游戏库理论软件

//
// $Source: /home/gambit/CVS/gambit/sources/poly/poly.imp,v $
// $Date: 2002/08/27 17:29:48 $
// $Revision: 1.3 $
//
// DESCRIPTION:
// Implementation of polynomial class
//
// 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 "poly.h"

//--------------------------------------------------------------------------
//                      class: polynomial<T>
//--------------------------------------------------------------------------

//--------------------------------------------------------------------------
//                 constructors and a destructor
//--------------------------------------------------------------------------

template <class T> polynomial<T>::polynomial(const polynomial<T>& x) 
: coeflist(x.coeflist)
{ 
}

template <class T> polynomial<T>::polynomial(const T& coeff, const int& deg) 
{ 
  if (coeff != (T)0) {
    for (int i = 0; i < deg; i++)
      coeflist += (T)0;
    coeflist += coeff;
  }
}

template <class T> 
polynomial<T>::polynomial(const gList<T>& coefficientlist) 
: coeflist(coefficientlist)
{ 
}

template <class T> 
polynomial<T>::polynomial(const gVector<T>& coefficientvector) 
: coeflist()
{ 
  for (int i = 1; i <= coefficientvector.Length(); i++)
    coeflist += coefficientvector[i];
}

template <class T> polynomial<T>::polynomial(const int deg) 
: coeflist()
{ 
  if (deg >= 0) { 
    //gout << "Error is polynomial int constructor.\n"; 
    exit(1);
  }
}

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

//--------------------------------------------------------------------------
//                             operators
//--------------------------------------------------------------------------

template <class T> polynomial<T>& 
                   polynomial<T>::operator= (const polynomial<T>& y)
{
  if (this!=&y) coeflist = y.coeflist;

  return *this;
}

template <class T> 
bool polynomial<T>::operator== (const polynomial<T>& 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 <class T> 
bool polynomial<T>::operator!= (const polynomial<T>& y) const
{
  return !(*this == y); 
}

template <class T> const T& polynomial<T>::operator[](const int index) const
{ 
  return coeflist[index+1];
}

template <class T> polynomial<T> 
                   polynomial<T>::operator+(const polynomial<T>& y) const
{
  if      (  Degree() < 0) return polynomial<T>(y); 
  else if (y.Degree() < 0) return polynomial<T>(*this);

  int max_degree;

  if (Degree() > y.Degree()) {max_degree =   Degree();}
  else                       {max_degree = y.Degree();}

  polynomial<T> sum;
  for (int i = 0; i <= max_degree; i++) {
    sum.coeflist += (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.Length() >= 1) && 
          (sum.coeflist[sum.coeflist.Length()] == (T)0) )
    sum.coeflist.Remove(sum.coeflist.Length());

  return sum;
}

template <class T> polynomial<T> 
                   polynomial<T>::operator-(const polynomial<T>& y) const
{
  return polynomial<T>(*this+(-y));
}

template <class T> polynomial<T> 
                   polynomial<T>::operator*(const polynomial<T> &y) const
{
  if      (Degree() == -1) return polynomial<T>(*this); 
  else if (y.Degree() == -1) return polynomial<T>(y);

  int tot_degree = Degree() + y.Degree();

  polynomial<T> product;
  for (int t = 0; t <= tot_degree; t++) 
    product.coeflist += (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 <class T> polynomial<T> 
                   polynomial<T>::operator/(const polynomial<T> &q) const
{
  assert(q.Degree() >= 0);

  polynomial<T> ans;
  polynomial<T> r = *this;
  while (r.Degree() >= q.Degree()) {
    polynomial<T> x(r.LeadingCoefficient()/q.LeadingCoefficient(),
		    r.Degree() - q.Degree());
    ans += x;
    r-=q*x;
  }
  return polynomial<T>(ans);
}

template <class T> polynomial<T>&
                   polynomial<T>::operator+=(const polynomial<T>& y) 
{
  return ((*this)=(*this)+y);
}

template <class T> polynomial<T>&
                   polynomial<T>::operator-=(const polynomial<T>& y)
{
  return ((*this)=(*this)-y);
}

template <class T> polynomial<T>&
                   polynomial<T>::operator*=(const polynomial<T>& y)
{
  return ((*this)=(*this)*y);
}

template <class T> polynomial<T>&
                   polynomial<T>::operator/=(const polynomial<T>& y)
{
  return ((*this)=(*this)/y);
}

template <class T> polynomial<T> 
                   polynomial<T>::operator%(const polynomial<T> &q) const
{
  assert (q.Degree() != -1);

  polynomial<T> ans;
  polynomial<T> r = *this;

  while (r.Degree() >= q.Degree())
    {
      polynomial<T> x(r.LeadingCoefficient()/q.LeadingCoefficient(),
		      r.Degree() - q.Degree());
      ans += x;
      r-=q*x;
    }
  return polynomial<T>(r);
}

template <class T> polynomial<T> polynomial<T>::operator- () const
{
  polynomial<T> negation;
  for (int i = 0; i <= Degree(); i++)
    negation.coeflist += -coeflist[i+1];

  return negation;
}

template <class T> polynomial<T> polynomial<T>::Derivative() const 
{
  if (Degree() <= 0) return polynomial<T>(-1);

  polynomial<T> derivative; 
  for (int i = 1; i <= Degree(); i++)
    derivative.coeflist += (T)i * coeflist[i+1];

  return derivative;
}

//--------------------------------------------------------------------------
//                             manipulation
//--------------------------------------------------------------------------

template <class T>  void polynomial<T>::ToMonic() 
{
  assert (!IsZero());

  T lc = LeadingCoefficient();
  for (int i = 1; i <= coeflist.Length(); i++)
    coeflist[i] /= lc;
}

template<class T>  
polynomial<gDouble>  polynomial<T>::TogDouble() const
{
  gList<gDouble> newcoefs;
  for (int i = 1; i <= coeflist.Length(); i++) {
    newcoefs += (gDouble)coeflist[i];
  }
  return polynomial<gDouble>(newcoefs);
}

//--------------------------------------------------------------------------
//                              information
//--------------------------------------------------------------------------

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

template <class T>  T polynomial<T>::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 <class T> int polynomial<T>::Degree() const
{
  return coeflist.Length() - 1;
}

template <class T> T polynomial<T>::LeadingCoefficient() const
{
  if (Degree() < 0) return (T)0;
  else              return coeflist[Degree() + 1];
}

template <class T> gList<T> polynomial<T>::CoefficientList() const
{
  return coeflist;
}

template <class T> 
polynomial<T> polynomial<T>::GcdWith(const polynomial<T>& that) const
{
  assert( !this->IsZero() && !that.IsZero() );

  polynomial<T> numerator(*this);
  numerator.ToMonic();
  polynomial<T> denominator(that);
  denominator.ToMonic();
  polynomial<T> remainder(numerator % denominator);

  while (!remainder.IsZero()) {
    remainder.ToMonic();
    numerator = denominator; 
    denominator = remainder; 
    remainder = numerator % denominator;
  }

  return denominator;
}

template <class T> 
bool polynomial<T>::IsQuadratfrei() const
{
  polynomial<T> Df(Derivative());
  if (Df.Degree() <= 0)             return true;
  if ( GcdWith(Df).Degree() <= 1 ) return true;
  else                              return false;
}

template <class T> 
bool polynomial<T>::CannotHaveRootsIn(const gInterval<T>& 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
  gList< polynomial<T> > derivs;
  derivs += Derivative();
  int i;
  for (i = 2; i <= Degree(); i++) derivs += 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,(long)max_index) < pow((T)base_of_max_index,(long)i))
	{ max_index = i; base_of_max_index = base; }
    }

  }
                             assert(base_of_max_index != (T)0);

  if ( (T)pow((T)I.Length(),(long)max_index) < (T)base_of_max_index ) 
    return true;
  else 
    return false;
}


template <class T> gList< gInterval<T> >
polynomial<T>::RootSubintervals(const gInterval<T>& I) const
{ assert ( Degree() >= 0 && IsQuadratfrei() );

  polynomial<T> Df = Derivative();

  gList< gInterval<T> > answer;

  gList< gInterval<T> > to_be_processed; to_be_processed += I;
  while (to_be_processed.Length() > 0) {

    gInterval<T> in_process = to_be_processed.Remove(to_be_processed.Length());

    if      ( EvaluationAt(in_process.LowerBound()) == (T)0 ) {

      if (Df.CannotHaveRootsIn(in_process)) {
	answer += in_process;
      }
      else {	
	to_be_processed += in_process.RightHalf();
	to_be_processed += in_process.LeftHalf();
      }
    }

    else if ( EvaluationAt(in_process.UpperBound()) == (T)0 ) {
      if (Df.CannotHaveRootsIn(in_process))  {
	if (in_process.UpperBound() == I.UpperBound()) {
	  answer += in_process;
	}
      }
      else {
	to_be_processed += in_process.RightHalf();
	to_be_processed += 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 += in_process;
	}
      }
      else {
	to_be_processed += in_process.RightHalf();
	to_be_processed += in_process.LeftHalf();
      }
    }

  }

  return answer;
}

template <class T> gInterval<T>
polynomial<T>::NeighborhoodOfRoot(const gInterval<T>& I, T& error) const
{