poly.imp

Gambit 是一个游戏库理论软件

  gList< gInterval<T> > intrvls; intrvls += I;
  while (intrvls[intrvls.Length()].Length() >= error) {
    if      (EvaluationAt(intrvls[intrvls.Length()].LowerBound()) == (T)0) 
      intrvls += gInterval<T>(intrvls[intrvls.Length()].LowerBound(),
			      intrvls[intrvls.Length()].LowerBound());

    else if (EvaluationAt(intrvls[intrvls.Length()].UpperBound()) == (T)0) 
      intrvls += gInterval<T>(intrvls[intrvls.Length()].UpperBound(),
			      intrvls[intrvls.Length()].UpperBound());

    else if ( (EvaluationAt(intrvls[intrvls.Length()].LowerBound()) >= (T)0 &&
	       EvaluationAt(intrvls[intrvls.Length()].Midpoint  ()) <= (T)0)   ||
              (EvaluationAt(intrvls[intrvls.Length()].LowerBound()) <= (T)0 &&
	       EvaluationAt(intrvls[intrvls.Length()].Midpoint  ()) >= (T)0) )
         intrvls += intrvls[intrvls.Length()].LeftHalf(); 

    else intrvls += intrvls[intrvls.Length()].RightHalf(); 
  }

  return intrvls[intrvls.Length()];

  // It is, perhaps, possible to speed this up, at least for double's
  // by introducing Newton's method.
}

template <class T> gList< gInterval<T> >
polynomial<T>::PreciseRootIntervals(const gInterval<T>& I, T& error) const
{
  gList< gInterval<T> > coarse = RootSubintervals(I);
  gList< gInterval<T> > fine;

  for (int i = 1; i <= coarse.Length(); i++) 
    fine += NeighborhoodOfRoot(coarse[i],error);

  return fine;
}

template <class T> gList<T>
polynomial<T>::PreciseRoots(const gInterval<T>& I, T& error) const
{
  gList<T> roots;

  polynomial<T> p(*this), factor(*this);

  while ( p.Degree() > 0 ) {

    int depth = 1;
    polynomial<T> probe(p.Derivative());
    polynomial<T> 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;
    gList< gInterval<T> > fine = factor.PreciseRootIntervals(I,error);
    for (int j = 1; j <= fine.Length(); 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 += approx;
    }
    factor = p;
  }

  return roots;
}

//--------------------------------------------------------------------------
//                                printing
//--------------------------------------------------------------------------

template <class T> void polynomial<T>::Output(gOutput &output) const
{
  if (Degree() < 0) output << "0";
  else {
    for (int t = 0; t <= Degree(); t++)
      if (coeflist[t+1] != (T)0) { 
	output << coeflist[t+1]; 
	if (t >= 1) output << "x";
	if (t >= 2) output << "^" << t;
	if (t < Degree()) output << " + "; }
  }
}

template <class T> gOutput &operator<<(gOutput &p_file,
				       const polynomial<T> &p_poly)
{
  p_poly.Output(p_file);
  return p_file;
}

//--------------------------------------------------------------------------
//                      class: complexpoly
//--------------------------------------------------------------------------

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

complexpoly::complexpoly(const complexpoly& x) 
: coeflist(x.coeflist)
{ 
}

complexpoly::complexpoly(const gComplex& coeff, const int& deg) 
{ 
  if (coeff != (gComplex)0) {
    for (int i = 0; i < deg; i++)
      coeflist += (gComplex)0;
    coeflist += coeff;
  }
}

 
complexpoly::complexpoly(const gList<gComplex>& coefficientlist) 
: coeflist(coefficientlist)
{ 
}

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

complexpoly::~complexpoly()
{
}

//--------------------------------------------------------------------------
//                             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 complexpoly(y); 
  else if (y.Degree() < 0) return complexpoly(*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 += (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.Length() >= 1) && 
          (sum.coeflist[sum.coeflist.Length()] == (gComplex)0) )
    sum.coeflist.Remove(sum.coeflist.Length());

  return sum;
}

complexpoly complexpoly::operator-(const complexpoly& y) const
{
  return complexpoly(*this+(-y));
}

complexpoly complexpoly::operator*(const complexpoly &y) const
{
  if      (Degree() == -1) return complexpoly(*this); 
  else if (y.Degree() == -1) return complexpoly(y);

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

  complexpoly product;
  for (int t = 0; t <= tot_degree; t++) 
    product.coeflist += (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 complexpoly(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 complexpoly(r);
}

complexpoly complexpoly::operator- () const
{
  complexpoly negation;
  for (int i = 0; i <= Degree(); i++)
    negation.coeflist += -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 += (gComplex)i * coeflist[i+1];

  return derivative;
}

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

void complexpoly::ToMonic() 
{
  assert (!IsZero());

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


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

bool complexpoly::IsZero() const
{
  if (coeflist.Length() == 0) return true;
  else                        return false;
}

gComplex complexpoly::EvaluationAt(const gComplex& arg) const
{
  gComplex answer;

  for (int i = 0; i <= Degree(); i++) {
    gComplex monom_val = (gComplex)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.Length() - 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;
}


gList<gComplex> complexpoly::Roots() const 
{
  assert (!IsZero());

  gList<gComplex> 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 += guess;

  gList<gComplex> lin_form_coeffs;
  lin_form_coeffs += guess;
  lin_form_coeffs += gComplex(-1.0,0.0);
  complexpoly linear_form(lin_form_coeffs);
  complexpoly quotient = *this/linear_form;
  answer += quotient.Roots();

  for (int i = 1; i <= answer.Length(); 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;
}

//--------------------------------------------------------------------------
//                                printing
//--------------------------------------------------------------------------

gOutput& operator << (gOutput& output, const complexpoly& x)
{
  if (x.Degree() < 0) output << "0";
  else {
    for (int t = 0; t <= x.Degree(); t++)
      if (x.coeflist[t+1] != (gComplex)0) { 
	output << x.coeflist[t+1]; 
	if (t >= 1) output << "x";
	if (t >= 2) output << "^" << t;
	if (t < x.Degree()) output << " + "; }
  }
  return output;
}