quiksolv.imp

Gambit 是一个游戏库理论软件

template <class T>
bool QuikSolv<T>::NewtonRootIsOnlyInRct(const gRectangle<gDouble>& r,
					gVector<gDouble>& point) const
{
  assert (NoEquations == System.Dmnsn());

  if ( NewtonRootInRectangle(r,point) ) {
    gSquareMatrix<gDouble> Df = TreesOfPartials.SquareDerivativeMatrix(point);
    if ( HasNoOtherRootsIn(r,point,Df.Inverse()) ) return true;
    else                                           return false;
  }
  else                                             return false;
}


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

template<class T> 
QuikSolv<T>& QuikSolv<T>::operator=(const QuikSolv<T> & rhs)
{
  assert (System == rhs.System);

  if (*this != rhs) {
    HasBeenSolved = rhs.HasBeenSolved;
    Roots         = rhs.Roots;
  }
  return *this;
}

template<class T>  
bool QuikSolv<T>::operator==(const QuikSolv<T> & rhs) const
{
    if (System        != rhs.System        ||
	HasBeenSolved != rhs.HasBeenSolved ||
	Roots         != rhs.Roots)
         return false;
    else return true;
}

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

//-------------------------------------------
//          Improve Accuracy of Root
//-------------------------------------------

template <class T> gVector<gDouble> 
QuikSolv<T>::NewtonPolishOnce(const gVector<gDouble>& point) const
{
  gVector<gDouble> oldevals = TreesOfPartials.ValuesOfRootPolys(point,
								NoEquations);
  gMatrix<gDouble> Df = TreesOfPartials.DerivativeMatrix(point,NoEquations);
  gSquareMatrix<gDouble> M(Df * Df.Transpose());
  
  gVector<gDouble> Del = - (Df.Transpose() * M.Inverse()) * oldevals;

  return point + Del;
}

template <class T> gVector<gDouble> 
QuikSolv<T>::SlowNewtonPolishOnce(const gVector<gDouble>& point) const
{
  gVector<gDouble> oldevals = TreesOfPartials.ValuesOfRootPolys(point,
								NoEquations);
  gMatrix<gDouble> Df = TreesOfPartials.DerivativeMatrix(point,NoEquations);
  gSquareMatrix<gDouble> M(Df * Df.Transpose());
  
  gVector<gDouble> Del = - (Df.Transpose() * M.Inverse()) * oldevals;

  bool done = false;
  while (!done) {
    gVector<gDouble> 
      newevals(TreesOfPartials.ValuesOfRootPolys(point + Del,NoEquations));
    if (newevals * newevals <= oldevals * oldevals) done = true;
    else for (int i = 1; i <= Del.Length(); i++) Del[i] /= 2;
  }

  return point + Del;
}

template<class T>   gVector<gDouble> 
QuikSolv<T>::NewtonPolishedRoot(const gVector<gDouble> &initial) const
{
  gList<gInterval<gDouble> > interval_list;
  for (int i = 1; i <= Dmnsn(); i++) 
    interval_list += gInterval<gDouble>(initial[i] - (gDouble)1, 
					initial[i] + (gDouble)1);
  gRectangle<gDouble> box(interval_list);
  gVector<gDouble> point(initial);
  if (!NewtonRootInRectangle(box,point))
    throw NewtonError();
  point = SlowNewtonPolishOnce(point);
  point = SlowNewtonPolishOnce(point);
  return point;
}

//-------------------------------------------
//          Check for Singularity
//-------------------------------------------

template<class T> bool
QuikSolv<T>::MightHaveSingularRoots() const
{
  assert (NoEquations == System.Dmnsn());

  gPoly<gDouble> newpoly = 
    gDoubleSystem.ReductionOf(  gDoubleSystem.DetOfDerivativeMatrix(),
			      *(gDoubleSystem.TermOrder()));

  if (newpoly.IsZero()) return true;

  gList<gPoly<gDouble> > newlist(gDoubleSystem.UnderlyingList());

  newlist += newpoly;
  gPolyList<gDouble> larger_system(gDoubleSystem.AmbientSpace(),
			     gDoubleSystem.TermOrder(),
			     newlist);
  gIdeal<gDouble> test_ideal(gDoubleSystem.TermOrder(),larger_system);

  return !(test_ideal.IsEntireRing());
}

//-------------------------------------------
//           The Central Calculation
//-------------------------------------------

template<class T> bool
QuikSolv<T>::FindRoots(const gRectangle<T>& r, const int& max_iterations) 
{
  assert (NoEquations == System.Dmnsn());

  int zero = 0;
  return FindCertainNumberOfRoots(r,max_iterations,zero);
}

template<class T> bool
QuikSolv<T>::FindCertainNumberOfRoots(const gRectangle<T>& r, 
				      const int& max_iterations,
				      const int& max_no_roots) 
{
  assert (NoEquations == System.Dmnsn());

  gList<gVector<gDouble> >* rootlistptr = new gList<gVector<gDouble> >();

  if (NoEquations == 0) {
    gVector<gDouble> answer(0);
    *rootlistptr += answer; 
    Roots = *rootlistptr;  
    HasBeenSolved = true;
    return true;
  }

  /* - Commmented out 7/11/00 because g3.nfg does not work
  gList<gVector<gDouble> > answer;
  bool done = PelicanRoots(r, answer);
  if (done) {
    Roots = answer;
    HasBeenSolved = true; 
    return true;     
  }
  */

  gArray<int> precedence(System.Length());  
            // Orders search for nonvanishing poly
  for (int i = 1; i <= System.Length(); i++) precedence[i] = i;

  int iterations = 0;

  int* no_found = new int(0);
  FindRootsRecursion(rootlistptr,
		     TogDouble(r), 
		     max_iterations, 
		     precedence, 
		     iterations,
		     max_no_roots,
		     no_found);

  if (iterations < max_iterations) { 
    Roots = *rootlistptr; 
    HasBeenSolved = true; 
    return true; 
  }

  return false;

/* - This is code that was once used to call the Grobner basis solver
     It will not work with T = gDouble or double, due to numerical instability
  gSolver<T> bigsolver(System.TermOrder(), System);

  if ( !bigsolver.IsZeroDimensional() ) return false;
  else {
    gList<gVector<gDouble> > rootlist;
    rootlist = bigsolver.Roots();
    gList<gVector<gDouble> > roots;
    for (int j = 1; j <= rootlist.Length(); j++)
      if (TogDouble(r).Contains(rootlist[j])) roots += rootlist[j];
    Roots = roots;
  }
  return true;
*/

}


template <class T> void // gList<gVector<gDouble> >  
QuikSolv<T>::FindRootsRecursion(      gList<gVector<gDouble> >* rootlistptr,
				const gRectangle<gDouble>& r, 
				const int& max_iterations,
				      gArray<int>& precedence,
				      int& iterations,
				const int& max_no_roots,
				      int* roots_found)    const
{
  assert (NoEquations == System.Dmnsn());

  // Check for user interrupt
  //  m_status.SetProgress(50.0);
  m_status.Get();

  if ( SystemHasNoRootsIn(r, precedence) ) 
    return;

  gVector<gDouble> point = r.Center();

  if ( NewtonRootIsOnlyInRct(r, point) ) {
    int i;
    for (i = NoEquations + 1; i <= System.Length(); i++)
      if (TreesOfPartials[i].ValueOfRootPoly(point) < (gDouble)0)
	return;

    bool already_found = false;
    for (i = 1; i <= rootlistptr->Length(); i++)
      if (point == (*rootlistptr)[i])
	already_found = true;
    if (!already_found) {
      *rootlistptr += point;
      (*roots_found)++;
    }
    return;
  }

  int N = r.NumberOfCellsInSubdivision();
  for (int i = 1; i <= N; i++)
    if (max_no_roots == 0 || *roots_found < max_no_roots) {
      if (iterations >= max_iterations) return;
      else {
	iterations++;
	FindRootsRecursion(rootlistptr,
			   r.SubdivisionCell(i),
			   max_iterations, 
			   precedence, 
			   iterations,
			   max_no_roots,
			   roots_found);
      } 
    }
  return; 
}



template <class T> const bool
QuikSolv<T>::ARootExistsRecursion(const gRectangle<gDouble>& r, 
					gVector<gDouble>& sample,
				  const gRectangle<gDouble>& smallrect, 
				        gArray<int>& precedence)    const
{
  if (smallrect.MaximalSideLength() == (gDouble)0.0) {
    sample = smallrect.Center();
    return true;
  }

  if ( SystemHasNoRootsIn(smallrect, precedence) ) 
    return false;		        

  gVector<gDouble> point = smallrect.Center();
  if (NewtonRootNearRectangle(smallrect,point))
    if (r.Contains(point)) {
      bool satisfies_inequalities(true);
      for (int i = NoEquations + 1; i <= System.Length(); i++)
	if (satisfies_inequalities)
	  if (TreesOfPartials[i].ValueOfRootPoly(point) < (gDouble)0)
	    satisfies_inequalities = false;
      if (satisfies_inequalities) {
	sample = point;
	return true;
      }
    }

  int N = smallrect.NumberOfCellsInSubdivision();
  for (int i = 1; i <= N; i++) 
    if (ARootExistsRecursion(r,
			     sample,
			     smallrect.SubdivisionCell(i),
			     precedence))
      return true;

  return false; 
}

template <class T> bool
QuikSolv<T>::ARootExists (const gRectangle<T>& r, 
			        gVector<gDouble>& sample) const
{
  if (NoEquations == 0) {
    gVector<gDouble> answer(0);
    sample = answer;
    return true;
  }

  gRectangle<gDouble> r_double(TogDouble(r));
  gArray<int> precedence(System.Length());  
            // Orders search for nonvanishing poly
  for (int i = 1; i <= System.Length(); i++) precedence[i] = i;

  return ARootExistsRecursion(r_double, sample, r_double, precedence);
}

//----------------------------------
//           Printing
//----------------------------------

template <class T> void QuikSolv<T>::Output(gOutput &output) const
{
  output << "The system is\n" << System;
  if (HasBeenSolved) 
    output << "It has been solved with roots:\n" << Roots;
  else output << "It has not been solved.";
}

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

template <class T> QuikSolv<T>::NewtonError::~NewtonError()
{ }

template <class T> gText QuikSolv<T>::NewtonError::Description(void) const 
{
  return "Newton method failed to polish approximate root";
}