//
// $Source: /cvsroot/gambit/gambit/sources/tools/enumpoly/gsolver.imp,v $
// $Date: 2006/08/17 03:05:00 $
// $Revision: 1.6 $
//
// DESCRIPTION:
// Implementation of class gSolver
//
// 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 "gsolver.h"

//---------------------------------------------------------------
//                      gSolver
//---------------------------------------------------------------

//---------------------------
// Constructor / Destructor
//---------------------------


template<class T>  gSolver<T>::gSolver(const term_order* Order,
				       const gPolyList<T>& Inputs)
: InputList(Inputs), TheIdeal(Order, Inputs)
{
}


template<class T>  gSolver<T>::gSolver(const gSolver<T>& given)
: InputList(given.InputList), TheIdeal(given.TheIdeal)
{
}

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

//---------------------------
//        Utilities
//---------------------------

template <class T> Gambit::List<gPoly<double> > gSolver<T>::BasisTogDouble() const
{
  Gambit::List<gPoly<double> > newlist;
  gPolyList<T> oldlist = TheIdeal.CanonicalBasis();
  for (int i = 1; i <= oldlist.Length(); i++) {
    
    newlist.Append(TogDouble(oldlist[i])); 
  }
  return newlist;
}

template <class T> Gambit::List<Gambit::Vector<double> >  
gSolver<T>::ContinuationSolutions(const Gambit::List<gPoly<double> >& list,
				  const int dmnsn,
				  const int curvar,
				  const Gambit::Vector<double>& knownvals)
{
  Gambit::List<Gambit::Vector<double> > answer;

  if (!list[1].IsZero()) {
    polynomial<double> rootpoly = list[1].UnivariateEquivalent(curvar);

    gInterval<double> root_interval((double)-10,(double)10); //  Ouch!!
    double error = (double) 0.0001;                           // Ditto!!

//DEBUG
//gout << "About to enter PreciseRoots ... \n";

    Gambit::List<double> rootlist = rootpoly.PreciseRoots(root_interval,error);

//DEBUG
//gout << "Just left PreciseRoots ... \n";


    for (int i = 1; i <= rootlist.Length(); i++) {
      Gambit::Vector<double> newvec(knownvals);
      newvec[curvar] = rootlist[i];
      if (curvar == dmnsn) answer.Append(newvec);
      else {
	Gambit::List<gPoly<double> > newlist;
	for (int j = 2; j <= list.Length(); j++) 
	  newlist.Append(list[j].EvaluateOneVar(curvar,rootlist[i]));
	answer += ContinuationSolutions(newlist,dmnsn,curvar+1,newvec);
      }
    }
  }
  
  else {   // (list[1].IsZero())
    Gambit::List<gPoly<double> > newlist;
    for (int j = 2; j <= list.Length(); j++) 
      newlist.Append(list[j]);
    answer += ContinuationSolutions(newlist,dmnsn,curvar,knownvals);
  }

  return answer;
}

//---------------------------
//        Information
//---------------------------

template<class T> bool gSolver<T>::IsZeroDimensional()
{
  return TheIdeal.ZeroDimensional();
}

/* - Old Implementation of Roots
template<class T> Gambit::List<Gambit::Vector<double> > gSolver<T>::Roots()
{
  Gambit::List<gPoly<double> >   dlist = BasisTogDouble();

  int dmnsn = InputList.AmbientSpace()->Dmnsn();

  Gambit::Vector<double> origin(dmnsn);
  for (int i = 1; i <= dmnsn; i++) origin[i] = (double)0; // Needless!

  Gambit::List<Gambit::Vector<double> > rootlist = ContinuationSolutions(dlist,
							    dmnsn,
							    1,
							    origin);
  return rootlist;
}
*/

template<class T> Gambit::List<Gambit::Vector<double> > gSolver<T>::Roots()
{
  Gambit::List<exp_vect> mon_bss = TheIdeal.MonomialBasis();

  Gambit::List<polynomial<T> > UnivariateRootEquations;
  int i;
  for (i = 1; i <= TheIdeal.Dmnsn(); i++) {
    Gambit::Matrix<T> CoefficientMatrix(mon_bss.Length()+1, mon_bss.Length());
    int j = 0; 
    bool done = false;
    while (!done) {
      exp_vect x_i_to_j(InputList.AmbientSpace(),i,j);
      gPoly<T> power(InputList.AmbientSpace(),x_i_to_j,(T)1,TheIdeal.Order());
      gPoly<T> reduced_power = TheIdeal.CanonicalBasis().
	                                ReductionOf(power,*(TheIdeal.Order()));
      for (int k = 1; k <= mon_bss.Length(); k++)
	CoefficientMatrix(j+1,k) = reduced_power.GetCoef(mon_bss[k]);

//DEBUG
//gout << "After updating, j = " << j
//  << "and CoefficientMatrix is\n" << CoefficientMatrix << "\n";

      Gambit::Matrix<T> SubMatrix(j+1,mon_bss.Length());
      for (int m = 1; m <= j+1; m++)
	for (int n = 1; n <= mon_bss.Length(); n++)
	  SubMatrix(m,n) = CoefficientMatrix(m,n);

//DEBUG
//gout << "Before entering Linear Combination, SubMatrix is\n" 
//  << SubMatrix << "\n";

      LinearCombination<T> Attempt(SubMatrix);
      if (Attempt.LastRowIsSpanned()) {
	polynomial<T> root_eqn_i(Attempt.LinearDependence());
	UnivariateRootEquations.Append(root_eqn_i);
	done = true;
      }
      j++;
    }
  }

  Gambit::List<polynomial<double> > ConvertedUnivariateRootEquations;
  for (i = 1; i <= UnivariateRootEquations.Length(); i++)
    ConvertedUnivariateRootEquations.Append(UnivariateRootEquations[i].TogDouble());

  Gambit::List<Gambit::Vector<double> > original;
  Gambit::List<Gambit::Vector<double> > revised;
  Gambit::Vector<double> zero(TheIdeal.Dmnsn());
  original.Append(zero);

    gInterval<double> root_interval((double)-10,(double)10); //  Ouch!!
    double error = (double) 0.0001;                           // Ditto!!

  for (i = 1; i <= TheIdeal.Dmnsn(); i++) {
    Gambit::List<double> roots_of_eqn_i =
      ConvertedUnivariateRootEquations[i].PreciseRoots(root_interval,error);
    for (int j = 1; j <= original.Length(); j++)
      for (int k = 1; k <= roots_of_eqn_i.Length(); k++) {
      Gambit::Vector<double> new_vec = original[j];
      new_vec[i] = roots_of_eqn_i[k];
      revised.Append(new_vec);
    }
    original = revised;
    revised = Gambit::List<Gambit::Vector<double> >();
  }

  Gambit::List<Gambit::Vector<double> > finished;
  Gambit::List<gPoly<double> > gDoublePolys = InputList.ListTogDouble();
  gPolyList<double> InputListTogDouble(TheIdeal.TheSpace(),
					TheIdeal.Order(),
					gDoublePolys);
  for (i = 1; i <= original.Length(); i++)
    if (InputListTogDouble.IsRoot(original[i])) finished.Append(original[i]);
  return finished;
}