realroot: include/realroot/solver_uv

00001 /*******************************************************************
00002  *   This file is part of the source code of the realroot kernel.
00003  *   Author(s): J.P. Pavone GALAAD, INRIA
00004  *              B. Mourrain GALAAD, INRIA
00005  ********************************************************************/
00006 #ifndef realroot_solver_sleeve_hpp
00007 #define realroot_solver_sleeve_hpp
00008 /********************************************************************/
00060 //====================================================================
00061 # include <realroot/IEEE754.hpp>
00062 # include <realroot/texp_bool.hpp>
00063 # include <realroot/array.hpp>
00064 # include <realroot/solver.hpp>
00065 # include <realroot/solver_binary.hpp>
00066 # include <realroot/ring_bernstein_tensor.hpp>
00067 # include <realroot/ring_monomial_tensor.hpp>
00068 
00069 # define TMPL    template<class C, class V>
00070 # define SOLVER  solver<C, Sleeve<V> >
00071 //====================================================================
00072 namespace mmx {
00073   template<class V> struct Sleeve{};
00074 
00075   template<class C>
00076   struct sleeve_rep {
00077     C *up, *dw;
00078     unsigned m_sz;
00079 
00080     sleeve_rep(unsigned n): up(new C[n]), dw(new C[n]), m_sz(n) {};
00081     unsigned size() const { return m_sz; }
00082   };
00083 
00084   struct Monomials;
00085   struct Bernstein;
00086   namespace let {
00087     template<class C, class D> void
00088     assign(sleeve_rep<C>& f, const polynomial<D, with<Bernstein> >& p) {
00089       assert(f.size()==(unsigned)(degree(p)+1));
00090       { 
00091         numerics::rounding<C> rnd(numerics::rnd_up());      
00092         for(unsigned i=0;i<f.size();i++) f.up[i]=as<C>(p[i]);
00093       }
00094       { 
00095         numerics::rounding<C> rnd(numerics::rnd_dw());      
00096         for(unsigned i=0;i<f.size();i++) f.dw[i]=as<C>(p[i]);
00097       }
00098     }
00099     template<class C, class D> void
00100     assign(sleeve_rep<C>& f, const polynomial<D, with<MonomialTensor> >& p) {
00101       assert(f.size()==(unsigned)(degree(p)+1));
00102       typedef typename kernelof<D>::T::integer  integer;
00103       typedef typename kernelof<D>::T::rational rational;
00104 
00105       binomials<integer> bn;
00106       unsigned ip=0, in=0;
00107       int d=degree(p);
00108       for(unsigned i=0;i< p.size();i++) {
00109         if(p[i]>p[ip]) {
00110           ip = i;
00111         } 
00112         if(p[i]<p[in]) in=i;
00113       }
00114       rational mx = (p[ip]>-p[in]?as<rational>(p[ip]):as<rational>(-p[in]));
00115       rational c;
00116       {
00117         numerics::rounding<C> rnd( numerics::rnd_up() );
00118         for ( int i = 0; i <= d; i ++ ) {
00119           c = as<rational>(p[i])/(mx*bn(d,i));
00120           f.up[i]=as<C>(c);
00121         }
00122       }
00123       {
00124         numerics::rounding<C> rnd( numerics::rnd_dw() );
00125         for ( int i = 0; i <= d; i ++ ) {
00126           c = as<rational>(p[i])/(mx*bn(d,i));
00127           f.dw[i]=as<C>(c);
00128         }
00129       } 
00130       //      polynom< ring<rational,Bernstein> > q;
00131       //      let::assign(q,p);
00132       //  assign(f,q);
00133     }
00134 
00135   }
00136   //====================================================================  
00137   template<class K, class Num>  struct binary_convert {};
00138   template<class K>
00139   struct binary_convert<K,Isolate>: public K {
00140 
00141     typedef Seq<Interval<typename K::ieee> >  Solutions;
00142     static data_t data;
00143 
00144     template<class output, class real> static inline
00145     void get(output& sol, const real& first = 0, const real& last =1)
00146     {
00147       typedef typename output::value_type root_t;
00148       typedef root_t                      as_root;
00149       //      unsigned s = o.size();
00150       real a,b;
00151       //      o.resize(s+2*size());
00152       real scale = last-first;
00153       for ( unsigned i = 0; i < data.nb_sol(); i ++ )
00154         {
00155           data_t::convert(a,data.m_res[i].a,first, scale);
00156           data_t::convert(b,data.m_res[i].b,first, scale);
00157           sol << as<root_t>(Interval<real>(a,b));
00158         };
00159     };
00160 
00161     static void set_precision(unsigned p) {};
00162 
00163     template<class output, class real> static inline
00164     void get(output & sol, const homography<real>& H)
00165     {
00166       typedef typename output::value_type root_t;
00167       typedef root_t                      as_root;
00168       typedef typename root_t::value_type bound_t;
00169       //      unsigned s = o.size();
00170       bound_t u,v;
00171       real dt= H.a*H.d-H.b*H.c;
00172       if(dt>0)
00173         for (unsigned i = 0; i < data.nb_sol(); i ++ )
00174           {
00175             data.convert(u,data.m_res[i].a,H);
00176             data.convert(v,data.m_res[i].b,H);
00177             sol << as_root(u,v);
00178           }
00179       else
00180         for ( int i = data.nb_sol()-1; i >=0 ; i --)
00181           {
00182             data.convert(u,data.m_res[i].b,H);
00183             data.convert(v,data.m_res[i].a,H);
00184             sol << as_root(u,v);
00185           }
00186     }
00187   };
00188   template <class K>  data_t binary_convert<K,Isolate>::data(true);
00189   //====================================================================
00190   template<class K>
00191   struct binary_convert<K,Approximate> : public K 
00192   {
00193     typedef unsigned unsigned_t;
00194     typedef Seq<typename K::ieee>  Solutions;
00195     typedef typename K::integer integer;
00196     static data_t data;
00197     
00198     static inline 
00199     integer binomial(const integer& n,const integer& p){return K::binomial(n,p);}
00200 
00201     template<class C> static inline C as_root(const C& a, const C& b) {return (a+b)/2;}
00202 
00203     template<class output, class real> static inline
00204     void get(output& sol, const real& first = 0, const real& last =1)
00205     {
00206       typedef typename output::value_type root;
00207       real a,b;
00208       real scale= last-first;
00209       for (unsigned i = 0; i < data.nb_sol(); i ++ )
00210         {
00211           data_t::convert(a,data.m_res[i].a, first, scale);
00212           data_t::convert(b,data.m_res[i].b, first, scale);
00213           sol << as<root>(Interval<real>(a,b));
00214         }
00215     } 
00216 
00217     static void set_precision(unsigned p) {};
00218 
00219     template<class output, class real> static inline
00220     void get(output& sol, const homography<real>& H)
00221     {
00222       //  typedef typename output::value_type as_root;
00223       //root_t u,v;
00224       unsigned s=data.nb_sol();
00225       std::vector<real> u(s);
00226       real v;
00227       real dt=H.a*H.d-H.b*H.c;
00228 
00229       if(dt>0)
00230         {
00231           {     
00232             numerics::rounding<real> rnd(numerics::rnd_dw());
00233             for ( unsigned i = 0; i < s; i ++ )
00234               data.convert(u[i],data.m_res[i].a,H);
00235           }
00236           {
00237             numerics::rounding<real> rnd(numerics::rnd_up());
00238             for ( unsigned i = 0; i < s; i ++ ) {
00239               data.convert(v,data.m_res[i].b,H);
00240               sol << as_root(u[i],v);
00241             }
00242           }
00243         }
00244       else
00245         {
00246           {     
00247             numerics::rounding<real> rnd(numerics::rnd_dw());
00248             for (int i = s-1; i >=0; i -- )
00249               data.convert(u[i],data.m_res[i].a,H);
00250           }
00251           {
00252             numerics::rounding<real> rnd(numerics::rnd_up());
00253             for (int i = s-1; i >=0; i -- ) {
00254               data.convert(v,data.m_res[i].b,H);
00255               sol << as_root(v,u[i]);
00256             }
00257           }
00258         }
00259     }
00260 
00261   };
00262 
00263   template <class K>  data_t binary_convert<K,Approximate>::data(false);
00264 
00265   //====================================================================
00267   template < class K>
00268   struct binary_sleeve_subdivision : public K {
00269 
00270     typedef typename K::integer  integer;
00271     typedef typename K::rational rational;
00272     //    typedef typename K::floating floating;
00273     //    typedef typename K::ieee     ieee;
00274 
00275     typedef double                      creal_t;
00276     typedef unsigned                    sz_t; 
00277     typedef binary_subdivision<K>       Base_t;
00278     typedef typename Base_t::unsigned_t unsigned_t;
00279     typedef res_t                       Domain_t;
00280   
00281     //    sleeve(): Base_t(1e-6){}
00282     //    sleeve(creal_t e): Base_t(e){}
00283 
00284     static inline void alloc( sz_t s, sz_t deep  ) 
00285     {
00286       K::data.alloc(s,2*s,deep); 
00287     };
00288     
00289     static void barre( char c, unsigned n )
00290     {
00291       char  _bar[ n+1 ];
00292       std::fill(_bar,_bar+n, c);
00293       _bar[n] = 0;
00294       std::cout << _bar << std::endl;
00295     };
00296     
00297     static void writebounds( creal_t * pup, creal_t * pdw, unsigned s )
00298     {
00299       Base_t::print(pup,s); std::cout << std::endl;
00300       Base_t::print(pdw,s); std::cout << std::endl;
00301     };
00302    
00303     static inline 
00304     bool glue( sz_t cup, sz_t cdw, int d )
00305     {
00306       if ( K::data.bckb() == K::data.m_dmn[d] )
00307         {
00308           K::data.bckb() = K::data.m_dmn[d];
00309           K::data.bckb().m += 1;
00310           if ( cup ) K::data.m_cup = cup;
00311           if ( cdw ) K::data.m_cdw = cdw;
00312           return true;
00313         };
00314       return false;
00315     };
00316     
00317     static inline 
00318     void mstore( sz_t cup, sz_t cdw, int d  ) 
00319     {
00320       if ( K::data.m_res.size() == 0 || !glue(cup,cdw,d)  )
00321         {
00322           K::data.mstore(d);
00323           K::data.m_cup = cup;
00324           K::data.m_cdw = cdw;
00325           //std::cout << "Push ";
00326         };
00327       //      writebck();
00328     };
00329   
00330 
00331     static inline
00332     void dwsplit( creal_t * r, creal_t * l, sz_t s )
00333     {
00334       numerics::rounding<creal_t> rnd( numerics::rnd_dw() );
00335       Base_t::split(r,l,s);
00336     };
00337 
00338     static inline 
00339     void upsplit( creal_t * r, creal_t * l, sz_t s )
00340     {
00341       numerics::rounding<creal_t> rnd( numerics::rnd_up() );
00342       Base_t::split(r,l,s);
00343     };
00344     
00345     static void Loop( bool isole = true ) {
00346       //      std::cout << "LOOP\n";
00347       numerics::rounding<creal_t> rnd(numerics::rnd_dw());
00348       creal_t * pup, * pdw;
00349       int d;
00350       sz_t a,s,cup,cdw;
00351       //      unsigned_t v;
00352       s = K::data.m_s;
00353       pup = K::data.m_data;
00354       pdw = K::data.m_data + s;
00355       K::data.m_dmn[0].m = 0;
00356       K::data.m_dmn[0].e = 0;
00357       for( a = 0, d = 0; d >= 0; d--, a -= 2*s   ) 
00358         {
00359           cup = Base_t::sgncnt(pup+a,s);
00360           cdw = Base_t::sgncnt(pdw+a,s);
00361           if ( cup || cdw )
00362             {
00363               if ( isole && cup == cdw && cup == 1 ){ 
00364                 //std::cout << " (b) \n";
00365                 K::data.sstore(d);
00366                 continue; 
00367               };
00368               if ( K::data.m_dmn[d].e == K::data.m_limit-1 ) { 
00369                 // std::cout << " (c) "<<cup<<" "<<cdw<<"\n";
00370                 if ( cup == cdw && cup == 1 )
00371                   K::data.sstore(d);
00372                 else
00373                   mstore(cup,cdw,d); 
00374                 continue; 
00375               };
00376               dwsplit(pdw+a,pdw+a+2*s,s);    
00377               upsplit(pup+a,pup+a+2*s,s);
00378               Base_t::split(K::data.m_dmn[d],K::data.m_dmn[d+1]);
00379               d += 2;
00380               a += 4*s;
00381               continue;
00382             } else {
00383             //std::cout<<"  (a) ";
00384             sz_t k =0, sv=0;
00385             for (  k = 0; k < K::data.m_s; k ++ ) {
00386               if ( ( pup[a+k] > 0 ) != ( pdw[a+k] > 0 ) ) {
00387                 sv++;
00388                 if(sv>1) { 
00389                   mstore(cup,cdw,d); 
00390                   break; 
00391                 };
00392               }
00393             }
00394             //std::cout<<std::endl;
00395           };
00396         };
00397       // std::cout << "END LOOP\n";
00398     };
00399     
00400 
00401     template<class input> static 
00402     void run_loop( const input& in, const creal_t& eps, const texp::true_t&   ) {
00403       //      binary_subdivision<K> slv(eps);
00404       sz_t prec = numerics::bitprec(eps);
00405       alloc(in.size(),prec);
00406       std::copy(in.begin(),in.end(),K::data.m_data);
00407       std::copy(in.begin(),in.end(),K::data.m_data+K::data.m_s);
00408       Loop(K::data.isole);
00409     };
00410     
00411     const Domain_t& operator[]( int i ) const  { return K::data.m_res[i]; };
00412     Domain_t&       operator[]( int i ) { return K::data.m_res[i]; };
00413     
00414   
00415     template<class output, class input> static inline
00416     void solve_bernstein( output& out, const input& in) {
00417       run(in,K::data.eps);
00418       //      solutions<typename output::value_type>::get(this,out);
00419     }
00420     template<class VECT,class POL,class Q> static inline 
00421     void init_pol(VECT& ubp, VECT& dbp, const POL& r,  unsigned sz, const Q& u, const Q&v);
00422 
00423 
00424     template<class input> static inline
00425     void run( const input& in, const creal_t& eps)
00426     { 
00427       typedef typename numerics::is_rounded<creal_t>::result_t round_t;
00428       run_loop(in,eps,round_t()); 
00429     }
00430     
00431     template<class C> static void 
00432     run( const sleeve_rep<C>& p) {
00433       sz_t prec = numerics::bitprec(K::data.eps);
00434       alloc(p.size(),prec);
00435       std::copy(p.up,p.up+p.size(),K::data.m_data);
00436       std::copy(p.dw,p.dw+p.size(),K::data.m_data+K::data.m_s);
00437       Loop(K::data.isole);
00438     }
00439 
00440     template<class input> static 
00441     void run( const input& up, const input& dw) {
00442       assert(up.size()==dw.size());
00443       sz_t prec = numerics::bitprec(K::data.eps);
00444       alloc(up.size(),prec);
00445       std::copy(up.begin(),up.end(),K::data.m_data);
00446       std::copy(dw.begin(),dw.end(),K::data.m_data+K::data.m_s);
00447       Loop(K::data.isole);
00448     }
00449     
00455     template <typename output, typename POL, typename Q> static 
00456     void solve(output& sol, const POL& r, const Q& u, const Q& v);
00457     
00458     //--------------------------------------------------------------------
00459     template<class output, class input,class real,class MTH> static inline
00460     void solve_bernstein( output& sol, 
00461                           const input& up, const input& dw, 
00462                           const real& u, const real& v, const MTH& mth)
00463     { 
00464       K::data.isole=mth.isole;
00465       run(up,dw);
00466       get(sol,u,v);
00467       //      solutions<typename output::value_type>::get(this,out,u,v);
00468     }
00469   
00470     //--------------------------------------------------------------------
00476     template <typename output, typename POL ,typename real, typename MTH> static 
00477     void solve_bernstein(output& sol, const POL& r, const real& u,const real& v, const MTH& mth )
00478     {
00479       typedef typename output::value_type    root_t;    
00480       using let::assign;
00481       
00482       unsigned sz= r.size();
00483       std::vector<creal_t> ubp( sz ), dbp( sz );
00484       {
00485         numerics::rounding<creal_t> rnd( numerics::rnd_up() );
00486         for(unsigned i=0;i<sz;i++) ubp[i] =as<creal_t>(r[i]);
00487       }
00488       {
00489         numerics::rounding<creal_t> rnd( numerics::rnd_dw() );
00490         for(unsigned i=0;i<sz;i++) dbp[i] =as<creal_t>(r[i]);
00491       }
00492       K::data.isole=mth.isole; 
00493       run(ubp,dbp);
00494       get(sol,u,v);
00495       //      solutions<typename output::value_type>::get(this, sol, u,v);
00496     }
00497   };
00498   //====================================================================
00499   template<class K>
00500   template<class VECT,class POL,class Q>
00501   void binary_sleeve_subdivision<K>::init_pol(VECT& ubp, VECT& dbp, const POL& r, unsigned sz,const Q& u, const Q&v)
00502   {
00503     using let::assign;
00504   
00505 #if 1
00506     //    PRINT_DEBUG("Conversion using exact arithmetic");
00507     rational U,V;
00508     assign(U,u);
00509     assign(V,v);
00510     //    std::vector<interval> bp( sz );
00511     std::vector<rational> bp(sz);
00512     univariate::convertm2b(bp,r,sz,U,V);
00513     rational mx = array::max_abs(bp);
00514     array::div(bp,mx);
00515     // unsigned i;
00516     //    for (i = 0; i < sz && bp[i].upper()-bp[i].lower()<as<floating>(1.e-16); i ++ );
00517     //    if(i<sz) std::cout<<"refine: "<<i<<" " <<(bp[i].upper()-bp[i].lower())<<std::endl;
00518     //    double l; assign(l, rational(V-U));
00519     //    std::cout<<"\t init: "<<sign_variation(bp)<<" "<<l<<std::endl;
00520     {
00521       numerics::rounding<creal_t> rnd( numerics::rnd_up() );
00522       for ( unsigned i = 0; i < sz; i ++ ) ubp[i]=as<creal_t>(bp[i]);
00523     }
00524     {
00525       numerics::rounding<creal_t> rnd( numerics::rnd_dw() );
00526       for ( unsigned i = 0; i < sz; i ++ ) dbp[i]=as<creal_t>(bp[i]);
00527     };
00528    
00529 #else
00530     //    PRINT_DEBUG("Using rounding arithmetic");
00531     O ud, vd;
00532     {
00533       numerics::rounding<creal_t> rnd( numerics::rnd_dw() );
00534       ud=as<creal_t>(u);
00535     };
00536     {
00537       numerics::rounding<creal_t> rnd( numerics::rnd_up() );
00538       vd=as<creal_t>(v);
00539       BEZIER::monomial_to_bezier(ubp,r,r.size(),ud,vd);
00540     };
00541     {
00542       numerics::rounding<creal_t> rnd( numerics::rnd_dw() );
00543       BEZIER::monomial_to_bezier(dbp,r,r.size(),ud,vd);
00544     };
00545 
00546     for ( unsigned i= 0; i < ubp.size(); i ++ )
00547       if ( dbp[i] > ubp[i] ) { 
00548         std::swap(dbp[i], ubp[i]);
00549       };
00550 #endif
00551 
00552   }
00553   //------------------------------------------------------------------
00554    template <class K>
00555   template <typename output, typename POL ,typename Q>
00556   void binary_sleeve_subdivision<K>::solve(output& sol, const POL& r, const Q& u, const Q& v)
00557   {
00558     unsigned sz= r.size();
00559     std::vector<double> ubp(sz), dbp(sz);
00560     init_pol(ubp,dbp,r,sz,u,v);
00561     run(ubp,dbp,true);
00562     //    run(ubp,dbp,this->eps,mth.isolate);
00563     //    solutions<typename output::value_type>::get(this,sol,u,v);
00564   }
00565 
00566   //====================================================================
00567   TMPL
00568   struct solver<C, Sleeve<V> > {
00569   
00570     typedef typename kernelof<C>::T          K;
00571     typedef typename binary_convert<K,V>::Solutions         Solutions;
00572     // Seq<Interval<typename K::ieee> >                    Solutions;
00573     typedef binary_sleeve_subdivision<binary_convert<K,V> > base_t;
00574 
00575 
00576     template<class POL> static void
00577     solve(Solutions& sol, const POL& p);
00578 
00579     template<class POL,class T> static void
00580     solve(Solutions& sol, const POL& p, const T& u, const T& v);
00581 
00582     template<class T, class U, class VRT> static void
00583     solve(Solutions& sol, const polynomial<U,VRT>& pl, const polynomial<U,VRT>& pu, const T& u, const T& v);
00584 
00585     template<class U,class VRT> static void
00586     solve(Solutions& sol, const polynomial<U,VRT>& pl, const polynomial<U,VRT>& pu);
00587 
00588 
00589     static inline void 
00590     set_precision(unsigned p) { binary_convert<K,V>::set_precision(p); }
00591 
00592   };
00593 
00596   TMPL
00597   template<class POL> void
00598   SOLVER::solve(Solutions& sol, const POL& p) {
00599       typedef typename K::ieee real_t;
00600       sleeve_rep<real_t> f(degree(p)+1); 
00601       let::assign(f,p);
00602       real_t c;
00603       for(unsigned i=1;i< f.size();i+=2) {
00604         c=f.up[i];f.up[i]=-f.dw[i];f.dw[i]=-c;
00605       }
00606       base_t::run(f);
00607 
00608       homography<real_t> H(1,0,1,-1);
00609       base_t::data.root_bound = - bound_root(p,Cauchy<real_t>());
00610       base_t::get(sol,H);
00611 
00612       for(unsigned i=1;i< f.size();i+=2) {
00613         c=f.up[i];f.up[i]=-f.dw[i];f.dw[i]=-c;
00614       }
00615       base_t::run(f);
00616       base_t::data.root_bound = - base_t::data.root_bound;
00617       H.c=-1;H.d=1;
00618       base_t::get(sol,H);
00619     }
00620   //--------------------------------------------------------------------
00621   TMPL
00622   template<class T, class U, class VRT> void
00623   SOLVER::solve(Solutions& sol, 
00624                 const polynomial<U,VRT>& pl, const polynomial<U,VRT>& pu,
00625                 const T& u, const T& v) {
00626     base_t::run(pl,pu);
00627     base_t::get(sol,u,v);
00628   }
00629 
00630   //--------------------------------------------------------------------
00631   TMPL
00632   template<class U, class VRT> void
00633   SOLVER::solve(Solutions& sol, 
00634                 const polynomial<U,VRT>& pl, const polynomial<U,VRT>& pu) {
00635     base_t::run(pl,pu);
00636     base_t::get(sol,0,1);
00637   }
00638 
00639   //--------------------------------------------------------------------
00640   TMPL
00641   template<class POL, class T> void
00642   SOLVER::solve(Solutions& sol, const POL& p,
00643                 const T& u, const T& v) {
00644     typedef typename POL::Ring                        Ring;
00645     typedef typename SOLVER::base_t      base_t;
00646     typedef typename SOLVER::Solutions Solutions;
00647     unsigned sz= p.size();
00648     std::vector<double> ubp(sz), dbp(sz);
00649     base_t::init_pol(ubp,dbp,p,sz,u,v);
00650     base_t::run(ubp,dbp);
00651     base_t::get(sol,u,v);
00652   }
00653   //--------------------------------------------------------------------
00654   template<class POL, class M, class B>
00655   typename solver<B, Sleeve<M> >::Solutions 
00656   solve (const POL& p, const Sleeve<M>& slv, const B& u, const B& v) {
00657     typedef solver<B, Sleeve<M> >                       Solver;
00658     typedef typename Solver::Solutions                  Solutions;
00659     Solutions sol;
00660     Solver::solve(sol,p, u,v);
00661     return sol;
00662   }
00663   //--x------------------------------------------------------------------
00664   template<class POL, class M>
00665   typename solver<typename POL::Scalar, Sleeve<M> >::Solutions 
00666   solve (const POL& p, const Sleeve<M>& slv) {
00667     typedef typename POL::Scalar                             Scalar;
00668     typedef solver<Scalar, Sleeve<M> >                       Solver;
00669     typedef typename solver<Scalar, Sleeve<M> >::Solutions   Solutions;
00670     Solutions sol;
00671     Solver::solve(sol,p);
00672     return sol;
00673   }
00674 
00675 } //namespace mmx
00676 //====================================================================
00677 # undef TMPL
00678 # undef SOLVER
00679 //====================================================================
00680 #endif  //realroot_solver_sleeve_hpp
00681
include/realroot/solver_uv_sleeve.hpp
