synaps/usolve/Sebastio.h

/********************************************************************
 *   This file is part of the source code of SYNAPS kernel.
 *   Author(s): B. Mourrain, H. Prieto, GALAAD, INRIA         
 *   $Id: Sebastio.h,v 1.1 2005/07/11 10:03:51 jppavone Exp $
 ********************************************************************/
#ifndef SYNAPS_SOLVE_SEBASTIO_H
#define SYNAPS_SOLVE_SEBASTIO_H
/********************************************************************/
#include <synaps/init.h>
#include <synaps/linalg/VECTOR.h>
#include <synaps/upoly/UPOL.h>
#include <synaps/upoly/UPolyQuot.h>

__BEGIN_NAMESPACE_SYNAPS

//--------------------------------------------------------------------
#define Tolerance 0.000001
#define Max_Iter 40
//--------------------------------------------------------------------
/* 
  The Sebastiao e Silva method is an iterative process for numerical
  approximation of the roots of a complex polynomial P.
  It is based on the power method, but instead of computing the successive 
  powers $f^{n}$,  one rather computes the sequence of powers $f^{2^{n}}$.
  See \cite{SeS41}, \cite{Card95}, for more details.
*/
//--------------------------------------------------------------------
struct Sebastiao {
  int nbroot;           //number of roots.
  unsigned nbiter;      //number of iterations.
  double   epsilon;     //tolerance.

  Sebastiao():nbroot(1),nbiter(40),epsilon(Tolerance){}
  Sebastiao(const Sebastiao & m):
    nbroot(m.nbroot),nbiter(m.nbiter),epsilon(m.epsilon){}
  Sebastiao(int n): 
    nbroot(n),nbiter(40),epsilon(Tolerance) {}
  Sebastiao(int n, unsigned nit): 
    nbroot(n),nbiter(nit),epsilon(Tolerance) {}
  Sebastiao(int n, unsigned nit, double eps): 
    nbroot(n),nbiter(nit),epsilon(eps) {}
};
//----------------------------------------------------------------------
template <class R> inline 
VAL<UPolyQuot<R>*> StartingValue(UPolDse<R> p, 
                                 const typename R::value_type x0)
{
  typename R::size_type degree = p.degree();
  UPolyQuot<R>* res = new UPolyQuot<R>(degree, Alloc());
  (*res)[degree-1]= (typename R::value_type)(1);
  for (typename R::size_type i = degree-1; i >=1; --i) {
    (*res)[i-1] = (*res)[i]*x0;
  }
  return VAL<UPolyQuot<R>*>(res);
}
//----------------------------------------------------------------------
template<class R>
VAL<UPolyQuot<R>*> IterFrom(const UPolyQuot<R> & p0, Sebastiao mth)
{
  assert(p0 !=0);
  
  typename R::size_type k=0, dim=UPolyQuot<R>::dim();

  RealOf<typename R::value_type>::TYPE bound(mth.epsilon);

  UPolyQuot<R>* un = new UPolyQuot<R>, u;

  u = p0;
  (*un) = u^2;
  (*un)/=VECTOR::norm(*un,dim);
  u    -=(*un);

  while ((VECTOR::norm(u,dim) > bound) && (k++ <= mth.nbiter))
    {
      u=(*un);
      (*un)=u^2;
      (*un)/=VECTOR::norm(*un,dim);
      u-=(*un);
    }
  #ifdef VERBATIM
   cout << " Number of iterations (in Sebastiao's method) : "<< k << endl;
  #endif
  if(k==mth.nbiter)
    cout <<"Warning: In Sebastio's method, maximum of iterations "
         <<k<<" reached"<<endl;
  return VAL<UPolyQuot<R>*>(un);
}
//----------------------------------------------------------------------
template <class R> inline 
typename R::value_type solveOne(UPolDse<R> & p, 
                             const typename R::value_type & x0, Sebastiao mth)
{
  typedef typename R::value_type C;
  UPolyQuot<R>::poly=p.rep;
  UPolyQuot<R>  u0;
  typename R::size_type dim=UPolyQuot<R>::dim();

  u0 = StartingValue(p, x0);  
  u0 = IterFrom(u0,mth);
  return u0[dim-2]/u0[dim-1];

}
//----------------------------------------------------------------------
template <class R> inline 
typename R::value_type* solve(UPolDse<R> & p, 
                              const typename R::value_type & x0, Sebastiao mth)
{
  typedef typename R::value_type C;
  typedef typename R::size_type  N;
  
  UPolyQuot<R>::poly = p.rep;
  UPolyQuot<R> sum(0),u0,u,un;
  N dim=UPolyQuot<R>::dim();
  
  int k;
  if(mth.nbroot<0) k=p.degree(); else k=mth.nbroot;

  C* x= new C[k];
  u0 = StartingValue(p, x0);

  for (int i = 0; i <k; i++)
    {
      if (i==0) 
        u=u0;
      else
        {u=u0; un=u*sum; u-=un;}

      un = IterFrom(u,mth);  
      x[i] = un[dim-2]/un[dim-1];

      /* First derivative of Poly estimated in the last root found,
         We use the coefficients of E, 
         which are the different powers of this root. */
      C deriv = 0;
      for (typename R::size_type j=1; j<=dim;j++)
        deriv+=UPolyQuot<R>::poly[j]*un[dim-j]*C(j);
      
      un/=deriv;
      sum += un;
    }
  return x;
}
//----------------------------------------------------------------------

__END_NAMESPACE_SYNAPS

#endif // SYNAPS_SOLVE_SEBASTIO_H