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