|
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): 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 ********************************************************************/ 00010 #ifndef realroot_solve_contfrac_hpp 00011 #define realroot_solve_contfrac_hpp 00012 00065 //==================================================================== 00066 #include <realroot/GMP.hpp> 00067 #include <realroot/univariate_bounds.hpp> 00068 #include <realroot/sign_variation.hpp> 00069 #include <realroot/Seq.hpp> 00070 #include <realroot/solver.hpp> 00071 #include <realroot/contfrac_intervaldata.hpp> 00072 #include <realroot/contfrac_lowerbound.hpp> 00073 00074 # define TMPL template<class C, class R, class V> 00075 # define SOLVER solver<C, ContFrac<R,V> > 00076 //==================================================================== 00077 namespace mmx { 00078 00079 template<class M> struct ContFrac{}; 00080 00081 //-------------------------------------------------------------------- 00082 template<class K, class B> struct continued_fraction_isolate: public K 00083 { 00084 typedef typename K::integer integer; 00085 typedef typename K::rational rational; 00086 typedef polynomial<integer, with<MonomialTensor> > Polynomial; 00087 typedef typename K::interval_rational root_t; 00088 typedef Seq<root_t> sol_t; 00089 00090 typedef IntervalData< integer, Polynomial> data; 00091 00092 static inline root_t as_root(const data& ID) 00093 { 00094 return as<root_t>(ID); 00095 } 00096 static inline bool is_empty(const data& ID) {return ID.s==0;} 00097 static inline bool is_good (const data& ID) {return ID.s==1;} 00098 static inline integer lower_bound (const Polynomial& p) {return B::lower_bound(p);} 00099 }; 00100 00101 00102 template<class K, class B> struct continued_fraction_approximate: public K 00103 { 00104 typedef typename K::integer integer; 00105 typedef typename K::rational rational; 00106 typedef typename K::interval_rational interval_rational; 00107 typedef polynomial<integer, with<MonomialTensor> > Polynomial; 00108 typedef interval_rational root_t; 00109 typedef Seq<root_t> sol_t; 00110 typedef IntervalData< integer, Polynomial> data; 00111 00112 static unsigned prec; 00113 static inline void set_precision(unsigned n){ prec = n;} 00114 static inline unsigned get_precision() { return prec; } 00115 00116 static inline root_t as_root(const data& ID) 00117 { 00118 return as<root_t>(ID); 00119 } 00120 00121 static inline bool is_empty(const data& ID) {return ID.s==0;} 00122 static inline bool is_good (const data& ID) 00123 { 00124 // std::cout << __FUNCTION__ <<" "<<ID.p<<" "<<ID.s<<std::endl; 00125 if(ID.s !=1) 00126 return false; 00127 else 00128 { 00129 integer N = ID.a*ID.d-ID.b*ID.c, D=ID.c*ID.d ; 00130 if(D==0) 00131 return false; 00132 if(N<0) N = -N; 00133 N<<=prec; 00134 if(D<0) D = -D; 00135 if(N<D) 00136 return true; 00137 else 00138 return false; 00139 } 00140 } 00141 static inline integer lower_bound(const Polynomial& p) {return B::lower_bound(p);} 00142 00143 }; 00144 00145 template<class K, class B> unsigned continued_fraction_approximate<K,B>::prec = 00146 solver_default_precision; 00147 00148 template < class K > 00149 struct continued_fraction_subdivision : public K { 00150 typedef typename K::integer RT; 00151 typedef typename K::integer integer; 00152 typedef typename K::rational FT; 00153 typedef typename K::rational rational; 00154 typedef typename K::floating floating; 00155 typedef typename K::root_t root_t; 00156 typedef typename K::Polynomial polynomial_integer; 00157 typedef polynomial< rational, with<MonomialTensor> > polynomial_rational; 00158 typedef polynomial< floating, with<MonomialTensor> > polynomial_floating; 00159 typedef typename K::data data; 00160 00161 static void set_precision(int prec) {K::set_precision(prec);} 00162 00163 static void 00164 solve_positive(Seq<root_t>& sol, const polynomial_integer& f, bool posneg) { 00165 // std::cout << __FUNCTION__ << " "<<f<<" "<<posneg<<std::endl; 00166 00167 typedef typename K::data IntervalData; 00168 00169 IntervalData ID(1, 0, 0, 1, f, 0); 00170 if (!posneg) 00171 { 00172 ID.a=-1; 00173 for (int i = 1; i <= degree(f); i+=2) ID.p[i] *= -1; 00174 } 00175 ID.s = sign_variation(ID.p); 00176 00177 if ( K::is_empty(ID) ) { return; } 00178 if ( K::is_good(ID) ) { 00179 //std::cout << "A) Sign variation: 1"<<std::endl;; 00180 FT B = bound_root( f, Cauchy<FT>()); 00181 if ( posneg ) 00182 sol << root_t( FT(0), B); 00183 else sol << root_t( FT(0), FT(-B)); 00184 return; 00185 } 00186 00187 solve_homography(sol,ID); 00188 } 00189 00190 template<class output> static void 00191 solve_homography(output& sol, const data& ID) { 00192 // std::cout << __FUNCTION__ << " "<<ID.p<<std::endl; 00193 00194 polynomial_integer X(1,0,1); 00195 //X[1] = RT(1); X[0]=RT(0); 00196 00197 int iters = 0; 00198 const RT One(1); 00199 FT Zero(0); 00200 // std::cout << "Polynomial is " << ID.p << std::endl; 00201 00202 typedef data IntervalData; 00203 std::stack< IntervalData > S; 00204 S.push( ID ); 00205 while ( !S.empty() ) { 00206 ++iters; 00207 IntervalData ID = S.top(); 00208 S.pop(); 00209 00210 // std::cout << "Polynomial is " << ID.p << std::endl; 00211 // Steps 3 - 4: Choose the method for computing the lower bound 00212 //-------------------------------------------------------------------- 00213 RT a = K::lower_bound(ID.p); 00214 // long k = Bd.lower_power_2( ID.p); 00215 // std::cout<< "\t Lower bound is "<< a<< std::endl; 00216 00217 if ( a > RT(1) ) { 00218 scale(ID,a); 00219 a = RT(1); 00220 } 00221 //-------------------------------------------------------------------- 00222 // Step 5 // 00223 if ( a >= RT(1) ) { 00224 shift(ID, a); 00225 // std::cout<<"Shift by "<<a<<std::endl; 00226 00227 if ( ID.p[0] == RT(0)) 00228 { 00229 div_x(ID.p,1); 00230 sol << root_t( to_FT( ID.b, ID.d), to_FT( ID.b, ID.d )); 00231 } 00232 ID.s = sign_variation( ID.p); 00233 } 00234 //int s= Kernel::init_shift(ID); 00235 //-------------------------------------------------------------------- 00236 if ( K::is_empty(ID) ) 00237 continue; 00238 else if( K::is_good(ID) ) { 00239 sol << K::as_root(ID); 00240 continue; 00241 } 00242 //-------------------------------------------------------------------- 00243 // Step 6 00244 IntervalData I1; 00245 shift_by_1( I1, ID); 00246 00247 unsigned long r = 0; 00248 00249 if (I1.p[0] == RT(0)) { 00250 // std::cout << "Zero at end point"<<std::endl;; 00251 sol << root_t( to_FT( I1.b, I1.d), to_FT(I1.b, I1.d) ); 00252 div_x(I1.p,1); 00253 r = 1; 00254 } 00255 I1.s = sign_variation( I1.p); 00256 00257 IntervalData I2; 00258 I2.s = ID.s - I1.s - r; 00259 reverse_shift_homography(I2,ID); 00260 00261 // std::cout <<I1.s <<" "<<I2.s<<std::endl; 00262 // Step 7; 00263 if ( !K::is_empty(I2) && !K::is_good(I2) ) { 00264 00265 //p2 = p2( 1 / (x+1)); 00266 reverse( I2.p, ID.p); 00267 shift_by_1(I2.p); 00268 00269 if ( I2.p[0] == 0) { 00270 div_x(I2.p,1); 00271 } 00272 I2.s = sign_variation( I2.p); 00273 } 00274 00275 // Step 8 00276 if ( I1.s < I2.s ) {std::swap( I1, I2); } 00277 00278 // Step 9 00279 if ( K::is_good(I1) ) { 00280 // std::cout << "C) Sign variation: 1"<<std::endl;; 00281 sol << K::as_root(I1); 00282 } else if ( !K::is_empty(I1) ) { 00283 S.push(I1); //IntervalData( I1.a, I1.b, I1.c, I1.d, I1.p, I1.s)); 00284 } 00285 // Step 10 00286 if ( K::is_good(I2) ) { 00287 // std::cout << "D) Sign variation: 1"<<std::endl;; 00288 sol << K::as_root(I2); 00289 } else if ( !K::is_empty(I2) ) { 00290 S.push(I2); //IntervalData( I2.a, I2.b, I2.c, I2.d, I2.p, I2.s)); 00291 } 00292 } // while 00293 // std::cout << "-------------------- iters: " << iters << std::endl; 00294 return; 00295 } 00296 00297 template < class output> inline static void 00298 solve(output& sol, const polynomial_integer& f, int mult=1) { 00299 // std::cout << __FUNCTION__ << " I "<<f<<std::endl; 00300 Seq < root_t > isolating_intervals; 00301 // Check if 0 is a root 00302 int idx; 00303 for (idx = 0; idx <= degree( f); ++idx) { 00304 if ( f[idx] != 0 ) break; 00305 } 00306 00307 polynomial_integer p; 00308 if ( idx == 0 ) { p = f; } 00309 else { 00310 p= polynomial_integer(1,degree(f) - idx); 00311 for (int j = idx; j <= degree(f); j++) 00312 p[j-idx] = f[j]; 00313 } 00314 00315 // std::cout << "Solving: " << p << std::endl; 00316 // Isolate the negative roots 00317 // for (int i = 1; i <= degree(p); i+=2) p[i] *= -1; 00318 solve_positive( isolating_intervals, p, false); 00319 // Is 0 a root? 00320 // std::cout << "ok negative" << std::endl; 00321 if (idx != 0) isolating_intervals << root_t( 0, 0); 00322 // for (int i = 1; i <= degree(p); i+=2) p[i] *= -1; 00323 solve_positive( isolating_intervals, p, true); 00324 // sort( isolating_intervals.begin(), isolating_intervals.end(), CompareFIT()); 00325 // std::cout << "ok positive" << std::endl; 00341 // std::cout << "now p: " << p << ", #nr: " << isolating_intervals.size() << std::endl; 00342 00343 // std::cout << "Done...." << std::endl; 00344 for (unsigned i = 0 ; i < isolating_intervals.size(); ++i) 00345 {sol << isolating_intervals[i];} 00346 // return sol; 00347 } 00348 00349 template < class output> inline static void 00350 solve(output& sol, const polynomial_rational& f) { 00351 // std::cout << __FUNCTION__ << " R "<<f<<std::endl; 00352 solve(sol, univariate::numer<polynomial_integer> (f)); 00353 } 00354 00355 template < class output> inline static void 00356 solve(output& sol, const polynomial_floating& f) { 00357 // std::cout << __FUNCTION__ << " R "<<f<<std::endl; 00358 polynomial_rational p(rational(0),degree(f)); 00359 00360 for(int i=0;i<degree(f)+1;i++) { 00361 p[i] = as<rational>(f[i]); 00362 } 00363 solve(sol,p); 00364 } 00365 00366 }; 00367 00368 //-------------------------------------------------------------------- 00369 template<class C> 00370 struct solver<C,ContFrac<Approximate> > { 00371 typedef typename texp::kernelof<C>::T K; 00372 typedef continued_fraction_subdivision< 00373 continued_fraction_approximate<K, 00374 AkritasBound<typename K::integer> > > 00375 base_t; 00376 typedef typename K::interval_rational interval_t; 00377 typedef Seq<interval_t> Solutions; 00378 typedef typename base_t::Polynomial polynomial_integer; 00379 typedef typename base_t::polynomial_rational polynomial_rational; 00380 typedef typename base_t::polynomial_floating polynomial_floating; 00381 00382 00383 static inline void set_precision(unsigned p) { 00384 continued_fraction_approximate<K,AkritasBound<typename K::integer> >::set_precision(p); 00385 } 00386 00387 template < class output > inline static void 00388 solve(output& sol, const polynomial_integer& f) { 00389 base_t::solve(sol, f); 00390 } 00391 00392 template < class output > inline 00393 static void solve(output& sol, const polynomial_rational& f) { 00394 base_t::solve(sol, f); 00395 } 00396 template < class output > inline static void 00397 solve(output& sol, const polynomial_floating& f) { 00398 base_t::solve(sol, f); 00399 } 00400 00401 template < class output > inline static void 00402 solve(output& sol, const polynomial_integer& f, int prec) { 00403 base_t::set_precision(prec); 00404 solve(sol, f); 00405 base_t::set_precision(solver_default_precision); 00406 } 00407 00408 template <class output> inline static void 00409 solve(output& sol, const polynomial_integer& f, 00410 const typename K::rational& u, const typename K::rational& v) { 00411 typedef typename K::integer integer; 00412 integer 00413 a= numerator(v), 00414 b= numerator(u), 00415 c= denominator(v), 00416 d= denominator(u); 00417 IntervalData<integer, polynomial_integer> D = 00418 as_interval_data(a,b,c,d,f); 00419 base_t::solve_homography(sol,D); 00420 } 00421 00422 template <class output> inline static void 00423 solve(output& sol, const polynomial_rational& f, 00424 const typename K::rational& u, const typename K::rational& v) { 00425 solve(sol, univariate::numer<polynomial_integer>(f), u, v); 00426 } 00427 00428 00429 }; 00430 00431 //==================================================================== 00432 template<class C> 00433 struct solver<C,ContFrac<Isolate> > { 00434 typedef typename texp::kernelof<C>::T K; 00435 typedef Seq<typename K::interval_rational> Solutions; 00436 typedef continued_fraction_subdivision< 00437 continued_fraction_isolate<K, 00438 AkritasBound<typename K::integer> > > base_t; 00439 00440 template < class output, class POL> inline static void 00441 solve(output& sol, const POL& f) { 00442 base_t::solve(sol,f); 00443 } 00444 }; 00445 00446 //--------------------------------------------------------------------- 00447 template<class POL, class M> 00448 typename solver<typename POL::Scalar, ContFrac<M> >::Solutions 00449 solve (const POL& p, const ContFrac<M>& slv) { 00450 typedef typename POL::Scalar Scalar; 00451 typedef solver<Scalar, ContFrac<M> > Solver; 00452 typedef typename solver<Scalar, ContFrac<M> >::Solutions Solutions; 00453 Solutions sol; 00454 Solver::solve(sol,p); 00455 return sol; 00456 } 00457 00458 //--------------------------------------------------------------------- 00459 template<class POL, class M, class B> 00460 typename solver<typename POL::Scalar, ContFrac<M> >::Solutions 00461 solve (const POL& p, const ContFrac<M>& slv, const B& b1, const B& b2) { 00462 typedef solver<B, ContFrac<M> > Solver; 00463 typedef typename solver<B, ContFrac<M> >::Solutions Solutions; 00464 Solutions sol; 00465 Solver::solve(sol,p,b1,b2); 00466 return sol; 00467 } 00468 //-------------------------------------------------------------------- 00469 00470 } //namespace mmx 00471 //==================================================================== 00472 # undef TMPL 00473 # undef SOLVER 00474 //==================================================================== 00475 #endif // _realroot_solve_contfrac_hpp_
