realroot: include/realroot/solver_uv_continued

00001 /*******************************************************************
00002  *   This file is part of the source code of the realroot kernel.
00003  *   Author(s):
00004  *        E. Tsigaridas <et@di.uoa.gr>
00005  *           Department of Informatics and Telecommunications
00006  *           University of Athens (Greece).
00007  *   Modification:
00008  *        B. Mourrain, GALAAD, INRIA
00009  *        P. Dobrowolski, Warsaw University of Technology
00010  ********************************************************************/
00011 #ifndef realroot_solve_contfrac_hpp
00012 #define realroot_solve_contfrac_hpp
00013 
00064 //====================================================================
00065 #include <realroot/GMP.hpp>
00066 #include <realroot/univariate_bounds.hpp>
00067 #include <realroot/sign_variation.hpp>
00068 #include <realroot/Seq.hpp>
00069 #include <realroot/solver.hpp>
00070 #include <realroot/contfrac_intervaldata.hpp>
00071 #include <realroot/contfrac_lowerbound.hpp>
00072 
00073 # define TMPL    template<class C, class R, class V>
00074 # define SOLVER  solver<C, ContFrac<R,V> >
00075 //====================================================================
00076 namespace mmx {
00077 
00078   template<class M> struct ContFrac{};
00079 
00080   //--------------------------------------------------------------------
00081   template<class K, class B, class Extended> struct continued_fraction_isolate;
00082 
00083   template<class K, class B> struct continued_fraction_isolate<K, B, false_t> : public K
00084   {
00085     typedef typename K::integer                                 integer;
00086     typedef typename K::rational                                rational;
00087     typedef polynomial<integer, with<MonomialTensor> >          Polynomial;
00088     typedef typename K::interval_rational                       root_t;
00089     typedef Seq<root_t>                                         sol_t;
00090     typedef IntervalData<integer, Polynomial>                   data;
00091 
00092     static inline root_t    as_root(const data& ID)             {return as<root_t>(ID);}
00093     static inline bool      is_empty(const data& ID)            {return ID.s==0;}
00094     static inline bool      is_good (const data& ID)            {return ID.s==1;}
00095     static inline integer   lower_bound(const Polynomial& p)    {return B::lower_bound(p);}
00096   };
00097 
00098   template<class K, class B> struct continued_fraction_isolate<K, B, true_t> : public K
00099   {
00100     typedef typename K::extended_integer                        integer;
00101     typedef typename K::extended_rational                       rational;
00102     typedef polynomial<integer, with<MonomialTensor> >          Polynomial;
00103     typedef typename K::interval_extended_rational              root_t;
00104     typedef Seq<root_t>                                         sol_t;
00105     typedef IntervalData<integer, Polynomial>                   data;
00106 
00107     static inline root_t    as_root(const data& ID)             {return as<root_t>(ID);}
00108     static inline bool      is_empty(const data& ID)            {return ID.s==0;}
00109     static inline bool      is_good (const data& ID)            {return ID.s==1;}
00110     static inline integer   lower_bound (const Polynomial& p)   {return B::lower_bound(p);}
00111   };
00112 
00113   template<class K, class B, class Extended> struct continued_fraction_approximate;
00114 
00115   template<class K, class B> struct continued_fraction_approximate<K, B, false_t> : public K
00116   {
00117     typedef typename K::integer                                 integer;
00118     typedef typename K::rational                                rational;
00119     typedef typename K::interval_rational                       interval_rational;
00120     typedef polynomial<integer, with<MonomialTensor> >          Polynomial;
00121     typedef interval_rational                                   root_t;
00122     typedef Seq<root_t>                                         sol_t;
00123     typedef IntervalData<integer, Polynomial>                   data;
00124 
00125     static unsigned prec;
00126     static inline void     set_precision(unsigned n){ prec = n;}
00127     static inline unsigned get_precision() { return prec; }
00128 
00129     static inline root_t as_root(const data& ID) { return as<root_t>(ID); }
00130 
00131     static inline bool is_empty(const data& ID) {return ID.s==0;}
00132     static inline bool is_good (const data& ID)
00133     {
00134       //      std::cout << __FUNCTION__ <<" "<<ID.p<<" "<<ID.s<<std::endl;
00135       if(ID.s !=1)
00136         return false;
00137       else
00138         {
00139           typename data::RT N = ID.a*ID.d-ID.b*ID.c, D=ID.c*ID.d ;
00140           if(D==0)
00141             return false;
00142           if(N<0) N = -N;
00143           N<<=prec;
00144           if(D<0) D = -D;
00145           if(N<D)
00146             return true;
00147           else
00148             return false;
00149         }
00150     }
00151 
00152     static inline integer           lower_bound(const Polynomial& p)            {return B::lower_bound(p);}
00153   };
00154 
00155   template<class K, class B> struct continued_fraction_approximate<K, B, true_t> : public K
00156   {
00157     typedef typename K::extended_integer                        integer;
00158     typedef typename K::extended_rational                       rational;
00159     typedef typename K::interval_extended_rational              interval_rational;
00160     typedef polynomial<integer, with<MonomialTensor> >          Polynomial;
00161     typedef interval_rational                                   root_t;
00162     typedef Seq<root_t>                                         sol_t;
00163     typedef IntervalData<integer, Polynomial>                   data;
00164 
00165     static unsigned prec;
00166     static inline void     set_precision(unsigned n){ prec = n;}
00167     static inline unsigned get_precision() { return prec; }
00168 
00169     static inline root_t as_root(const data& ID) { return as<root_t>(ID); }
00170 
00171     static inline bool is_empty(const data& ID) {return ID.s==0;}
00172     static inline bool is_good (const data& ID)
00173     {
00174       //      std::cout << __FUNCTION__ <<" "<<ID.p<<" "<<ID.s<<std::endl;
00175       if(ID.s !=1)
00176     return false;
00177       else
00178     {
00179           typename data::RT N = ID.a*ID.d-ID.b*ID.c, D=ID.c*ID.d ;
00180       if(D==0)
00181         return false;
00182       if(N<0) N = -N;
00183           N<<=prec;
00184       if(D<0) D = -D;
00185       if(N<D)
00186         return true;
00187       else
00188         return false;
00189     }
00190     }
00191 
00192     static inline integer           lower_bound(const Polynomial& p)            {return B::lower_bound(p);}
00193   };
00194 
00195   template<class K, class B> unsigned continued_fraction_approximate<K,B,false_t>::prec =
00196               solver_default_precision;
00197   template<class K, class B> unsigned continued_fraction_approximate<K,B,true_t>::prec =
00198         solver_default_precision;
00199 
00200   template < class K >
00201   struct continued_fraction_subdivision : public K {
00202     typedef typename K::integer       RT;
00203     typedef typename K::rational      FT;
00204     typedef typename K::integer       integer;
00205     typedef typename K::rational      rational;
00206     typedef typename K::floating      floating;
00207     typedef typename K::root_t        root_t;
00208     typedef typename K::Polynomial    polynomial_integer;
00209     typedef polynomial< rational, with<MonomialTensor> >  polynomial_rational;
00210     typedef polynomial< floating, with<MonomialTensor> >  polynomial_floating;
00211     typedef typename K::data          data;
00212 
00213     static void set_precision(int prec) {K::set_precision(prec);}
00214 
00215     template<class poly>
00216     static void
00217     solve_positive(Seq<root_t>& sol, const poly& f, bool posneg) {
00218       //      std::cout << __FUNCTION__ << " "<<f<<" "<<posneg<<std::endl;
00219 
00220       typedef typename K::data IntervalData;
00221 
00222             IntervalData ID(1, 0, 0, 1, f, 0);
00223             if (!posneg)
00224             {
00225                 ID.a=-1;
00226                 for (int i = 1; i <= degree(f); i+=2) ID.p[i] *= -1;
00227             }
00228             ID.s = sign_variation(ID.p);
00229 
00230             if ( K::is_empty(ID) ) { return; }
00231             if ( K::is_good(ID)  ) {
00232                 //std::cout << "A) Sign variation: 1"<<std::endl;;
00233                 FT B = bound_root( f, Cauchy<FT>());
00234                 if ( posneg )
00235                      sol << root_t( FT(0), B);
00236                 else sol << root_t( FT(0), FT(-B));
00237                 return;
00238             }
00239 
00240             solve_homography(sol,ID, poly());
00241     }
00242 
00243     template<class output, class data_type, class poly> static void
00244     solve_homography(output& sol, const data_type& ID, const poly &) {
00245         //            std::cout << __FUNCTION__ << " "<<ID.p<<std::endl;
00246 
00247 //        typename data_type::Poly X(1,0,1);
00248       //X[1] = RT(1); X[0]=RT(0);
00249 
00250       int iters  = 0;
00251       //const RT One(1);
00252       //FT Zero(0);
00253       //    std::cout << "Polynomial is " << ID.p << std::endl;
00254 
00255       typedef data_type IntervalData;
00256       std::stack< IntervalData > S;
00257       S.push( ID );
00258       while ( !S.empty() ) {
00259         ++iters;
00260         IntervalData ID = S.top();
00261         S.pop();
00262 
00263         //std::cout << "Polynomial is " << ID.p << std::endl;
00264         // Steps 3 - 4: Choose the method for computing the lower bound
00265         //--------------------------------------------------------------------
00266         RT a = K::lower_bound(ID.p);
00267         //    long k = Bd.lower_power_2( ID.p);
00268         //std::cout<< "\t Lower bound is "<< a<< std::endl;
00269 
00270         if ( a > RT(1) ) {
00271             scale(ID,a);
00272             a = RT(1);
00273         }
00274         //--------------------------------------------------------------------
00275         // Step 5 //
00276         if ( a >= RT(1) ) {
00277           shift(ID, a);
00278           //            std::cout<<"Shift by "<<a<<std::endl;
00279 
00280           if ( ID.p[0] == RT(0))
00281             {
00282               div_x(ID.p,1);
00283               sol << root_t( to_FT( ID.b, ID.d), to_FT( ID.b, ID.d ));
00284             }
00285           ID.s = sign_variation( ID.p);
00286         }
00287         //int s= Kernel::init_shift(ID);
00288         //--------------------------------------------------------------------
00289         if ( K::is_empty(ID) )
00290           continue;
00291         else if( K::is_good(ID) ) {
00292           sol << K::as_root(ID);
00293           continue;
00294         }
00295         //--------------------------------------------------------------------
00296         // Step 6
00297         IntervalData I1;
00298         shift_by_1( I1, ID);
00299 
00300         unsigned long r  = 0;
00301 
00302         if (I1.p[0] == RT(0)) {
00303           //    std::cout << "Zero at end point"<<std::endl;;
00304           sol << root_t( to_FT( I1.b, I1.d), to_FT(I1.b, I1.d) );
00305           div_x(I1.p,1);
00306           r = 1;
00307         }
00308         I1.s = sign_variation( I1.p);
00309 
00310         IntervalData I2;
00311         I2.s = ID.s - I1.s - r;
00312         reverse_shift_homography(I2,ID);
00313 
00314         // std::cout <<I1.s <<" "<<I2.s<<std::endl;
00315         // Step 7;
00316         if ( !K::is_empty(I2) && !K::is_good(I2) ) {
00317 
00318           //p2 = p2( 1 / (x+1));
00319           reverse( I2.p, ID.p);
00320           shift_by_1(I2.p);
00321 
00322           if ( I2.p[0] == 0) {
00323             div_x(I2.p,1);
00324           }
00325           I2.s = sign_variation( I2.p);
00326         }
00327 
00328         // Step 8
00329         if ( I1.s < I2.s ) {std::swap( I1, I2); }
00330 
00331         // Step 9
00332         if ( K::is_good(I1) ) {
00333           //    std::cout << "C) Sign variation: 1"<<std::endl;;
00334           sol << K::as_root(I1);
00335         } else if ( !K::is_empty(I1) ) {
00336           S.push(I1); //IntervalData( I1.a, I1.b, I1.c, I1.d, I1.p, I1.s));
00337         }
00338         // Step 10
00339         if ( K::is_good(I2) ) {
00340           //    std::cout << "D) Sign variation: 1"<<std::endl;;
00341           sol <<  K::as_root(I2);
00342         } else if ( !K::is_empty(I2) ) {
00343           S.push(I2); //IntervalData( I2.a, I2.b, I2.c, I2.d, I2.p, I2.s));
00344         }
00345       } // while
00346       //    std::cout << "-------------------- iters: " << iters << std::endl;
00347       return;
00348     }
00349 
00350     // solve polynomial_integer or polynomial_extended_integer
00351     template < class output, class poly> inline static void
00352     solve_polynomial(output& sol,  const poly& f, int  mult=1) {
00353       //      std::cout << __FUNCTION__ << " I "<<f<<std::endl;
00354       Seq < root_t > isolating_intervals;
00355       // Check if 0 is a root
00356       int idx;
00357       for (idx = 0; idx <= degree( f); ++idx) {
00358         if ( f[idx] != 0 ) break;
00359       }
00360 
00361       poly p;
00362       if ( idx == 0 ) { p = f; }
00363       else {
00364         p= poly(1,degree(f) - idx);
00365         for (int j = idx; j <= degree(f); j++)
00366           p[j-idx] = f[j];
00367       }
00368 
00369       // std::cout << "Solving: " << p << std::endl;
00370       // Isolate the negative roots
00371       //    for (int i = 1; i <= degree(p); i+=2) p[i] *= -1;
00372       solve_positive( isolating_intervals, p, false);
00373       // Is 0 a root?
00374       // std::cout << "ok negative" << std::endl;
00375       if (idx != 0) isolating_intervals << root_t( 0, 0);
00376       //    for (int i = 1; i <= degree(p); i+=2) p[i] *= -1;
00377       solve_positive( isolating_intervals, p, true);
00378       // sort( isolating_intervals.begin(), isolating_intervals.end(), CompareFIT());
00379       // std::cout << "ok positive" << std::endl;
00395       // std::cout << "now p: " << p << ",  #nr: " << isolating_intervals.size() << std::endl;
00396 
00397       //      std::cout << "Done...." << std::endl;
00398       for (unsigned i = 0 ; i < isolating_intervals.size(); ++i)
00399         {sol << isolating_intervals[i];}
00400       //    return sol;
00401     }
00402 
00403     template < class output> inline static void
00404     solve(output& sol,  const polynomial_integer& f) {
00405       solve_polynomial(sol, f);
00406     }
00407 
00408     template < class output> inline static void
00409     solve(output& sol,  const polynomial_rational& f) {
00410       //      std::cout << __FUNCTION__ << " R "<<f<<std::endl;
00411       solve_polynomial(sol, univariate::numer<polynomial_integer> (f));
00412     }
00413 
00414     template < class output> inline static void
00415     solve(output& sol,  const polynomial_floating& f) {
00416       //      std::cout << __FUNCTION__ << " R "<<f<<std::endl;
00417       polynomial_rational p(rational(0),degree(f));
00418 
00419       for(int i=0;i<degree(f)+1;i++) {
00420         p[i] = as<rational>(f[i]);
00421       }
00422       //solve_polynomial(sol,p);
00423       solve_polynomial(sol, univariate::numer<polynomial_integer> (p));
00424     }
00425 
00426   };
00427 
00428   //--------------------------------------------------------------------
00429   template<class C, class Extended> struct solver_approximate_traits;
00430 
00431   template<class C>
00432   struct solver_approximate_traits<C, false_t> {
00433     typedef typename texp::kernelof<C>::T K;
00434     typedef continued_fraction_subdivision<
00435               continued_fraction_approximate<K,
00436                 AkritasBound<typename K::integer>, false_t > >  base_t;
00437     typedef typename K::interval_rational                       interval_t;
00438   };
00439 
00440   template<class C>
00441   struct solver_approximate_traits<C, true_t> {
00442     typedef typename texp::kernelof<C>::T K;
00443     typedef continued_fraction_subdivision<
00444               continued_fraction_approximate<K,
00445                 AkritasBound<typename K::extended_integer>, true_t > >  base_t;
00446     typedef typename K::interval_extended_rational                      interval_t;
00447   };
00448 
00449   template<class C>
00450   struct solver<C,ContFrac<Approximate> > {
00451     typedef typename is_extended<C>::T              Extended;
00452     typedef solver_approximate_traits<C, Extended>  Traits;
00453     typedef typename Traits::K                      K;
00454     typedef typename Traits::base_t                 base_t;
00455     typedef typename Traits::interval_t             interval_t;
00456     typedef Seq<interval_t>                         Solutions;
00457     typedef typename base_t::Polynomial             polynomial_integer;
00458     typedef typename base_t::polynomial_rational    polynomial_rational;
00459     typedef typename base_t::polynomial_floating    polynomial_floating;
00460 
00461     static inline void set_precision(unsigned p) {
00462        continued_fraction_approximate<K,AkritasBound<typename K::integer>,Extended>::set_precision(p);
00463     }
00464 
00465     template < class output > inline static void
00466     solve(output& sol, const polynomial_integer& f) {
00467        base_t::solve(sol, f);
00468     }
00469 
00470     template < class output > inline
00471     static void solve(output& sol, const polynomial_rational& f) {
00472       base_t::solve(sol, f);
00473     }
00474 
00475     template < class output > inline static void
00476     solve(output& sol, const polynomial_floating& f) {
00477        base_t::solve(sol, f);
00478     }
00479 
00480     template < class output > inline static void
00481     solve(output& sol, const polynomial_integer& f, int prec) {
00482       base_t::set_precision(prec);
00483       solve(sol, f);
00484       base_t::set_precision(solver_default_precision);
00485     }
00486 
00487     template <class output> inline static void
00488     solve(output& sol, const polynomial_integer& f,
00489           const  typename K::rational& u, const  typename K::rational& v) {
00490       typedef typename K::integer integer;
00491       integer
00492         a= numerator(v),
00493         b= numerator(u),
00494         c= denominator(v),
00495         d= denominator(u);
00496       IntervalData<integer, polynomial_integer> D =
00497         as_interval_data(a,b,c,d,f);
00498       base_t::solve_homography(sol,D,polynomial_integer());
00499     }
00500 
00501     template <class output> inline static void
00502     solve(output& sol, const polynomial_rational& f,
00503           const  typename K::rational& u, const  typename K::rational& v) {
00504       solve(sol, univariate::numer<polynomial_integer>(f), u, v);
00505     }
00506 
00507 
00508   };
00509 
00510   //====================================================================
00511   template<class C, class Extended> struct solver_isolate_traits;
00512 
00513   template<class C>
00514   struct solver_isolate_traits<C, false_t> {
00515     typedef typename texp::kernelof<C>::T K;
00516     typedef continued_fraction_subdivision<
00517               continued_fraction_isolate<K,
00518                 AkritasBound<typename K::integer>, false_t > >  base_t;
00519     typedef typename K::interval_rational                       interval_t;
00520   };
00521 
00522   template<class C>
00523   struct solver_isolate_traits<C, true_t> {
00524     typedef typename texp::kernelof<C>::T K;
00525     typedef continued_fraction_subdivision<
00526               continued_fraction_isolate<K,
00527                 AkritasBound<typename K::extended_integer>, true_t > >  base_t;
00528     typedef typename K::interval_extended_rational                      interval_t;
00529   };
00530 
00531   template<class C>
00532   struct solver<C,ContFrac<Isolate> > {
00533     typedef typename is_extended<C>::T          Extended;
00534     typedef solver_isolate_traits<C,Extended>   Traits;
00535     typedef typename Traits::base_t             base_t;
00536     typedef typename Traits::interval_t         interval_t;
00537     typedef Seq<interval_t>                     Solutions;
00538 
00539     template < class output, class POL> inline static void
00540     solve(output& sol, const POL& f) {
00541       base_t::solve(sol,f);
00542     }
00543   };
00544 
00545  //---------------------------------------------------------------------
00546   template<class POL, class M>
00547   typename solver<typename POL::Scalar, ContFrac<M> >::Solutions
00548   solve (const POL& p, const ContFrac<M>& slv) {
00549     typedef typename POL::Scalar                             Scalar;
00550     typedef solver<Scalar, ContFrac<M> >                     Solver;
00551     typedef typename solver<Scalar, ContFrac<M> >::Solutions Solutions;
00552     Solutions sol;
00553     Solver::solve(sol,p);
00554     return sol;
00555   }
00556 
00557   //---------------------------------------------------------------------
00558   template<class POL, class M, class B>
00559   typename solver<typename POL::Scalar, ContFrac<M> >::Solutions
00560   solve (const POL& p, const ContFrac<M>& slv, const B& b1, const B& b2) {
00561     typedef solver<B, ContFrac<M> >                     Solver;
00562     typedef typename solver<B, ContFrac<M> >::Solutions Solutions;
00563     Solutions sol;
00564     Solver::solve(sol,p,b1,b2);
00565     return sol;
00566   }
00567   //--------------------------------------------------------------------
00568 
00569 } //namespace mmx
00570 //====================================================================
00571 # undef TMPL
00572 # undef SOLVER
00573 //====================================================================
00574 #endif // _realroot_solve_contfrac_hpp_
include/realroot/solver_uv_continued_fraction.hpp
