# synaps/usolve/Sebastio.h

```00001 /********************************************************************
00002  *   This file is part of the source code of SYNAPS kernel.
00003  *   Author(s): B. Mourrain, H. Prieto, GALAAD, INRIA
00004  *   \$Id: Sebastio.h,v 1.1 2005/07/11 10:03:51 jppavone Exp \$
00005  ********************************************************************/
00006 #ifndef SYNAPS_SOLVE_SEBASTIO_H
00007 #define SYNAPS_SOLVE_SEBASTIO_H
00008 /********************************************************************/
00009 #include <synaps/init.h>
00010 #include <synaps/linalg/VECTOR.h>
00011 #include <synaps/upoly/UPOL.h>
00012 #include <synaps/upoly/UPolyQuot.h>
00013
00014 __BEGIN_NAMESPACE_SYNAPS
00015
00016 //--------------------------------------------------------------------
00017 #define Tolerance 0.000001
00018 #define Max_Iter 40
00019 //--------------------------------------------------------------------
00020 /*
00021   The Sebastiao e Silva method is an iterative process for numerical
00022   approximation of the roots of a complex polynomial P.
00023   It is based on the power method, but instead of computing the successive
00024   powers \$f^{n}\$,  one rather computes the sequence of powers \$f^{2^{n}}\$.
00025   See \cite{SeS41}, \cite{Card95}, for more details.
00026 */
00027 //--------------------------------------------------------------------
00033 struct Sebastiao {
00034   int nbroot;           //number of roots.
00035   unsigned nbiter;      //number of iterations.
00036   double   epsilon;     //tolerance.
00037
00038   Sebastiao():nbroot(1),nbiter(40),epsilon(Tolerance){}
00039   Sebastiao(const Sebastiao & m):
00040     nbroot(m.nbroot),nbiter(m.nbiter),epsilon(m.epsilon){}
00041   Sebastiao(int n):
00042     nbroot(n),nbiter(40),epsilon(Tolerance) {}
00043   Sebastiao(int n, unsigned nit):
00044     nbroot(n),nbiter(nit),epsilon(Tolerance) {}
00045   Sebastiao(int n, unsigned nit, double eps):
00046     nbroot(n),nbiter(nit),epsilon(eps) {}
00047 };
00048 //----------------------------------------------------------------------
00052 template <class R> inline
00053 VAL<UPolyQuot<R>*> StartingValue(UPolDse<R> p,
00054                                  const typename R::value_type x0)
00055 {
00056   typename R::size_type degree = p.degree();
00057   UPolyQuot<R>* res = new UPolyQuot<R>(degree, Alloc());
00058   (*res)[degree-1]= (typename R::value_type)(1);
00059   for (typename R::size_type i = degree-1; i >=1; --i) {
00060     (*res)[i-1] = (*res)[i]*x0;
00061   }
00062   return VAL<UPolyQuot<R>*>(res);
00063 }
00064 //----------------------------------------------------------------------
00070 template<class R>
00071 VAL<UPolyQuot<R>*> IterFrom(const UPolyQuot<R> & p0, Sebastiao mth)
00072 {
00073   assert(p0 !=0);
00074
00075   typename R::size_type k=0, dim=UPolyQuot<R>::dim();
00076
00077   RealOf<typename R::value_type>::TYPE bound(mth.epsilon);
00078
00079   UPolyQuot<R>* un = new UPolyQuot<R>, u;
00080
00081   u = p0;
00082   (*un) = u^2;
00083   (*un)/=VECTOR::norm(*un,dim);
00084   u    -=(*un);
00085
00086   while ((VECTOR::norm(u,dim) > bound) && (k++ <= mth.nbiter))
00087     {
00088       u=(*un);
00089       (*un)=u^2;
00090       (*un)/=VECTOR::norm(*un,dim);
00091       u-=(*un);
00092     }
00093   #ifdef VERBATIM
00094    cout << " Number of iterations (in Sebastiao's method) : "<< k << endl;
00095   #endif
00096   if(k==mth.nbiter)
00097     cout <<"Warning: In Sebastio's method, maximum of iterations "
00098          <<k<<" reached"<<endl;
00099   return VAL<UPolyQuot<R>*>(un);
00100 }
00101 //----------------------------------------------------------------------
00108 template <class R> inline
00109 typename R::value_type solveOne(UPolDse<R> & p,
00110                              const typename R::value_type & x0, Sebastiao mth)
00111 {
00112   typedef typename R::value_type C;
00113   UPolyQuot<R>::poly=p.rep;
00114   UPolyQuot<R>  u0;
00115   typename R::size_type dim=UPolyQuot<R>::dim();
00116
00117   u0 = StartingValue(p, x0);
00118   u0 = IterFrom(u0,mth);
00119   return u0[dim-2]/u0[dim-1];
00120
00121 }
00122 //----------------------------------------------------------------------
00127 template <class R> inline
00128 typename R::value_type* solve(UPolDse<R> & p,
00129                               const typename R::value_type & x0, Sebastiao mth)
00130 {
00131   typedef typename R::value_type C;
00132   typedef typename R::size_type  N;
00133
00134   UPolyQuot<R>::poly = p.rep;
00135   UPolyQuot<R> sum(0),u0,u,un;
00136   N dim=UPolyQuot<R>::dim();
00137
00138   int k;
00139   if(mth.nbroot<0) k=p.degree(); else k=mth.nbroot;
00140
00141   C* x= new C[k];
00142   u0 = StartingValue(p, x0);
00143
00144   for (int i = 0; i <k; i++)
00145     {
00146       if (i==0)
00147         u=u0;
00148       else
00149         {u=u0; un=u*sum; u-=un;}
00150
00151       un = IterFrom(u,mth);
00152       x[i] = un[dim-2]/un[dim-1];
00153
00154       /* First derivative of Poly estimated in the last root found,
00155          We use the coefficients of E,
00156          which are the different powers of this root. */
00157       C deriv = 0;
00158       for (typename R::size_type j=1; j<=dim;j++)
00159         deriv+=UPolyQuot<R>::poly[j]*un[dim-j]*C(j);
00160
00161       un/=deriv;
00162       sum += un;
00163     }
00164   return x;
00165 }
00166 //----------------------------------------------------------------------
00167
00168 __END_NAMESPACE_SYNAPS
00169
00170 #endif // SYNAPS_SOLVE_SEBASTIO_H
00171
```

SYNAPS DOCUMENTATION