|
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 * V. Sharma, GALAAD, INRIA 00010 ********************************************************************/ 00011 00012 00013 #ifndef _realroot_SOLVE_CONTFRAC_HPP_ 00014 #define _realroot_SOLVE_CONTFRAC_HPP_ 00015 00016 #include <realroot/Interval.hpp> 00017 //#include <realroot/traits.hpp> 00018 #include <realroot/upoldse.hpp> 00019 #include <realroot/upoldse_bound.hpp> 00020 #include <realroot/sign.hpp> 00021 #include <realroot/Seq.hpp> 00022 00023 #include <realroot/contfrac_intervaldata.hpp> 00024 #include <realroot/contfrac_lowerbound.hpp> 00025 00026 00027 00028 namespace mmx { 00029 00030 00031 namespace meta { 00032 00033 class Kioustelidis_bound_1 00034 { 00035 00036 public: 00039 Kioustelidis_bound_1() {} 00040 00041 00042 template <typename NT> 00043 long upper( const upoldse<NT>& p) 00044 { 00045 return foo( p.begin(), p.end()); 00046 // return foo( begin( p), end( p)); 00047 } 00048 00049 template <typename NT> 00050 long lower( const upoldse<NT>& p) 00051 { 00052 return -foo( p.rbegin(), p.rend()); 00053 } 00054 00055 template <class IT> 00056 long foo( IT first, IT last) 00057 { 00058 typedef typename std::iterator_traits<IT>::difference_type size_type; 00059 00060 size_type deg = std::distance( first, last); 00061 assert( deg > 0 ); 00062 00063 int s = sign( *(last-1)); 00064 long lan = bitsize( *(last-1)); 00065 00066 bool first_time = true; 00067 long maxpow = 0; 00068 00069 int i = 1; 00070 for ( IT it = (last-2); it != first-1; --it, ++i) { 00071 if ( s * sign( *it) < 0) { 00072 long p = bitsize( abs( *it)) - lan - 1; 00073 long q = p / i; 00074 if ( first_time ) { 00075 first_time = false; 00076 maxpow = q + 2; 00077 } else { 00078 maxpow = std::max( maxpow, q + 2); 00079 } 00080 } 00081 } 00082 return maxpow; 00083 } 00084 }; 00085 00086 template <typename T> 00087 void times_2_to_k( T& r, const T& a, long k = 1) 00088 { 00089 r = a << k; 00090 00091 } 00092 00093 template <typename T> 00094 unsigned long int bitsize( const T& a) 00095 { 00096 return std::log( static_cast<double>(a)); 00097 } 00098 00099 #define __GMP_EXPR ::__gmp_expr<T, U> 00100 #define TMPL template <class T, class U> inline 00101 00102 TMPL 00103 mpz_class times_2_to_k( const __GMP_EXPR& a, long k = 1) 00104 { 00105 mpz_class b; 00106 mpz_mul_2exp( b.get_mpz_t(), a.get_mpz_t(), k ); 00107 return b; 00108 } 00109 00110 TMPL 00111 void times_2_to_k( __GMP_EXPR& r, const __GMP_EXPR& a, long k = 1) 00112 { 00113 mpz_mul_2exp( r.get_mpz_t(), a.get_mpz_t(), k ); 00114 } 00115 00116 TMPL 00117 unsigned long int bitsize( const __GMP_EXPR& a) 00118 { 00119 return mpz_sizeinbase( mpz_class( a).get_mpz_t(), 2 ); 00120 } 00121 00122 #undef TMPL 00123 #undef __GMP_EXPR 00124 00125 00126 00127 template < typename RT, typename Poly > inline 00128 void do_scale( IntervalData<RT,Poly>& ID, long k) 00129 { 00130 assert( k > 0 ); 00131 00132 for (int i = 0; i <= degree( ID.p); ++i) { 00133 times_2_to_k( ID.p[i], ID.p[i], k*i); 00134 } 00135 times_2_to_k( ID.a, ID.a, k); 00136 times_2_to_k( ID.c, ID.c, k); 00137 } 00138 00139 } //meta 00140 00141 00142 template<class C> typename meta::rationalof<C>::T to_rational(const C& a, const C& b) 00143 { 00144 typedef typename meta::rationalof<C>::T rational; 00145 return (rational(a)/b); 00146 } 00147 00149 template < class NT , class LowerBound = AkritasBound<NT> > 00150 struct ContFrac_t 00151 { 00152 typedef NT RT; 00153 typedef typename meta::rationalof<RT>::T FT; 00154 typedef Interval<FT> FIT; 00155 typedef Interval<FT> root_t; 00156 typedef upoldse<NT> Poly; 00157 typedef LowerBound bound; 00158 typedef ContFrac_t self_t; 00159 }; 00160 00161 00162 template < class NT, class LB > inline 00163 Seq< typename ContFrac<NT,LB>::root_t > 00164 solve( const typename ContFrac<NT>::Poly& f, ContFrac<NT,LB>) 00165 { 00166 typedef ContFrac<NT,LB> K; 00167 typedef typename K::root_t root_t; 00168 00169 Seq<root_t> sol; 00170 // std::cout << "**In this" << std::endl; 00171 solve_contfrac(f, std::back_inserter(sol.rep()), K()); 00172 // std::stable_sort( sol.begin(), sol.end(), ALGEBRAIC::Refine_compare()); 00173 00174 return sol; 00175 } 00176 00177 00178 template < typename K > 00179 void 00180 CF_positive( const typename K::Poly& f, Seq< typename K::FIT>& RL, bool posneg, K) 00181 { 00182 // std::cout << __FUNCTION__ << std::endl; 00183 typedef typename K::RT RT; 00184 typedef typename K::FT FT; 00185 typedef typename K::FT rational_t; 00186 typedef typename K::FIT FIT; 00187 typedef typename K::Poly Poly; 00188 00189 typedef IntervalData< RT, Poly> IntervalData; 00190 00191 meta ::Kioustelidis_bound_1 bound; 00192 00193 00194 int iters = 0; 00195 const RT One( 1); 00196 FT Zero(0); 00197 00198 // FT B = bound_root( f, Cauchy<RT>()); 00199 00200 FT B = meta::times_2_to_k( RT( 1), bound.upper( f)); 00201 00202 00203 unsigned long s = sign_variation(f); 00204 00205 if ( s == 0 ) { return; } 00206 if ( s == 1 ) { 00207 // std::cout << "A) Sign variation: 1"<<std::endl;; 00208 if ( posneg ) 00209 { 00210 // std::cout<<to_FIT( FT(0), B)<<std::endl;; 00211 RL.push_back( to_FIT( FT(0), B)); 00212 } 00213 else RL.push_back( to_FIT( FT(0), FT(-B))); 00214 return; 00215 } 00216 Poly X(2, AsSize()); X[1] = RT(1); X[0]=RT(0); 00217 00218 std::stack< IntervalData > S; 00219 if ( posneg) S.push( IntervalData( 1, 0, 0, 1, f, s)); 00220 else S.push( IntervalData( -1, 0, 0, 1, f, s)); 00221 00222 while ( !S.empty() ) { 00223 ++iters; 00224 IntervalData ID = S.top(); 00225 S.pop(); 00226 00227 //std::cout << "Polynomial is " << ID.p << std::endl; 00228 // Steps 3 - 4: Choose the method for computing the lower bound 00229 00230 00231 //typename K::bound Bd; 00232 // Cauchy<RT> Bd; 00233 //HongBound<RT> Bd; 00234 00235 00236 // RT M = Bd.upper(ID.p); 00237 // scale(ID,M); 00238 // reverse(ID); 00239 // translate_by_1(ID); 00240 // reverse(ID); 00241 00242 // Interval<rational_t> q=Bd.range(ID.p); 00243 // std::cout<<"Sign var: "<< sign_variation(ID.p)<<" "; 00244 // RT a, b; 00245 // if(q<0) 00246 // { 00247 // continue; 00248 // } 00249 // else if(q!=0) 00250 // { 00251 // a=numerator(q.lower()); 00252 // b=denominator(q.lower()); 00253 // inv_scale(ID,b); 00254 // shift(ID, a); 00255 // scale(ID,b); 00256 // let::assign(a,rational_t(q.upper()-q.lower())); 00257 // shift(ID,a); 00258 // a=One; 00259 // } 00260 00261 // std::cout<<"("<< sign_variation(ID.p)<<")"<<std::endl;; 00262 00263 00264 // RT a = Bd.lower(ID.p); 00265 //long k = bound.lower( ID.p); 00266 long k = AkritasBound<RT>().lower_power_2( ID.p); 00267 // std::cout << "k: " << k << std::endl; 00268 // std::cout<< "\t Lower bound is "<< a<< std::endl; 00269 00270 if (k > 0) { 00271 meta::do_scale( ID, k); 00272 } 00273 00274 // if ( a > One ) { scale(ID,a); a = One; } 00275 00276 // Step 5 // 00277 // if ( a >= One ) 00278 if ( k >= 0 ) { 00279 shift_by_1(ID); 00280 //shift(ID, a); 00281 00282 if ( ID.p[0] == RT(0)) { 00283 RL.push_back( to_FIT( to_rational( ID.b, ID.d), to_rational( ID.b, ID.d ))); 00284 ID.p /= X; 00285 } 00286 ID.s = sign_variation( ID.p); 00287 if ( ID.s == 0 ) { continue; } 00288 if ( ID.s == 1 ) { 00289 // std::cout << "B) Sign variation: 1"<<std::endl; 00290 RL.push_back( to_FIT<K>( ID.a, ID.b, ID.c, ID.d, B)); 00291 continue; 00292 } 00293 } 00294 00295 // Step 6 00296 IntervalData I1; 00297 shift_by_1( I1, ID); 00298 00299 unsigned long r = 0; 00300 00301 if (I1.p[0] == RT(0)) { 00302 // std::cout << "Zero at end point"<<std::endl;; 00303 RL.push_back( to_FIT( to_rational( I1.b, I1.d), to_rational(I1.b, I1.d) )); 00304 I1.p /= X; 00305 r = 1; 00306 } 00307 I1.s = sign_variation( I1.p); 00308 00309 IntervalData I2; 00310 00311 I2.s = ID.s - I1.s - r; 00312 00313 I2.a = ID.b; 00314 I2.b = ID.a + ID.b; 00315 I2.c = ID.d; 00316 I2.d = ID.c + ID.d; 00317 00318 // Step 7; 00319 if ( I2.s > 1 ) { 00320 //p2 = p2( 1 / (x+1)); 00321 reverse( I2.p, ID.p); 00322 shift_by_1( I2.p); 00323 00324 if ( I2.p[0] == 0) { 00325 I2.p /= X; 00326 } 00327 I2.s = sign_variation( I2.p); 00328 } 00329 00330 // Step 8 00331 if ( I1.s < I2.s ) { 00332 std::swap( I1, I2); 00333 } 00334 00335 // Step 9 00336 if ( I1.s == 1 ) { 00337 // std::cout << "C) Sign variation: 1"<<std::endl;; 00338 RL.push_back( to_FIT<K>( I1.a, I1.b, I1.c, I1.d, B)); 00339 } else if ( I1.s > 1) { 00340 S.push( IntervalData( I1.a, I1.b, I1.c, I1.d, I1.p, I1.s)); 00341 } 00342 00343 // Step 10 00344 if ( I2.s == 1 ) { 00345 // std::cout << "D) Sign variation: 1"<<std::endl;; 00346 RL.push_back( to_FIT<K>( I2.a, I2.b, I2.c, I2.d, B)); 00347 } else if ( I2.s > 1) { 00348 S.push( IntervalData( I2.a, I2.b, I2.c, I2.d, I2.p, I2.s)); 00349 } 00350 } // while 00351 // std::cout << "-------------------- iters: " << iters << std::endl; 00352 return; 00353 } 00354 00355 // A continued fraction algorithm that combines the Descartes method and 00356 // Akritas' algorithm. The basic idea is to take the better of performing 00357 // a Taylor shift or subdividing into half. 00358 template < typename K > 00359 void 00360 MCF_positive( const typename K::Poly& f, Seq< typename K::root_t>& RL, bool posneg, K) 00361 { 00362 // std::cout << __FUNCTION__ << std::endl; 00363 typedef typename K::RT RT; 00364 typedef typename K::FT FT; 00365 typedef typename K::root_t FIT; 00366 typedef typename K::Poly Poly; 00367 00368 typedef IntervalData< RT, Poly> IntervalData; 00369 00370 int iters = 0; 00371 //std::cout <<"Polynomial is "<< f << std::endl; 00372 00373 FT B = bound_root( f, Cauchy<FT>()); 00374 // std::cout << "Upper bound is "<< B << " Positive or negative " << posneg << std::endl; 00375 unsigned long s = sign_variation( f); 00376 if ( s == 0 ) { return; } 00377 if ( s == 1 ) { 00378 if ( posneg ) RL.push_back( to_FIT( FT(0), B)); 00379 else RL.push_back( to_FIT( FT(0), FT(-B))); 00380 return; 00381 } 00382 00383 Poly X(2, AsSize()); X[1] = RT(1); 00384 00385 std::stack< IntervalData > S; 00386 if ( posneg) S.push( IntervalData( 1, 0, 0, 1, f, s)); 00387 else S.push( IntervalData( -1, 0, 0, 1, f, s)); 00388 00389 IntervalData ID, I1, I2;// ID is the data associated with the parent and 00390 // I1, I2 are the data associated with its children 00391 00392 RT lowerBound; 00393 FT midpoint; // The point corresponding to the midpoint of the interval 00394 // with endpoints (a/c), (b/d). 00395 FT temp; // A temporary variable 00396 unsigned long mult = 0; // multiplicity of an exact root. 00397 00398 while ( !S.empty() ) { 00399 ++iters; 00400 ID = S.top(); 00401 S.pop(); 00402 // std::cout <<"Transformation is " << ID.a << ":" << ID.b << ":"<< ID.c << ":" << ID.d << std::endl; 00403 00404 // Steps 3 - 4 00405 lowerBound = CauchyMinimumBound( ID.p); 00406 std::cout << "lower bound is "<< lowerBound << std::endl; 00407 00408 if(ID.c == 0){// then ID.a is not equal to zero 00409 midpoint = FT(ID.d * B - ID.b*sign(ID.a), 2* abs(ID.a)); 00410 //std::cout << "Midpoint = " << midpoint << std::endl; 00411 }else if (ID.d > ID.c ){ 00412 midpoint = FT(ID.d, ID.c); 00413 }else{} 00414 00415 if(midpoint < 1) midpoint = 1; 00416 std::cout << "Midpoint = " << midpoint << std::endl; 00417 00418 if(lowerBound >= midpoint){ 00419 if ( lowerBound > 1 ) { 00420 UPOLDSE::scale( ID.p, ID.p, lowerBound); 00421 ID.a *= lowerBound; 00422 ID.c *= lowerBound; 00423 lowerBound = 1; 00424 } 00425 00426 if ( lowerBound == 1 ) { 00427 shift_by_1( ID.p, ID.p);// If we scaled above then this is required 00428 ID.b += ID.a; // to bring the origin to lowerBound 00429 ID.d += ID.c; 00430 00431 if ( ID.p[0] == RT(0)) { 00432 RL.push_back( to_FIT( FT( ID.b, ID.d), FT( ID.b, ID.d))); 00433 ID.p /= X; 00434 } 00435 ID.s = sign_variation( ID.p); 00436 if ( ID.s == 0 ) { continue; } 00437 if ( ID.s == 1 ) { 00438 RL.push_back( to_FIT<K>( ID.a, ID.b, ID.c, ID.d, B)); 00439 continue; 00440 } 00441 } 00442 00443 shift_by_1( I1.p, ID.p); 00444 I1.a = ID.a; 00445 I1.b = ID.a + ID.b; 00446 I1.c = ID.c; 00447 I1.d = ID.c + ID.d; 00448 00449 if (I1.p[0] == RT(0)) { 00450 RL.push_back( to_FIT( FT( I1.b, I1.d), FT( I1.b, I1.d))); 00451 I1.p /= X; 00452 mult = 1; 00453 } 00454 I1.s = sign_variation( I1.p); 00455 00456 I2.s = ID.s - I1.s - mult; 00457 I2.a = ID.b; 00458 I2.b = ID.a + ID.b; 00459 I2.c = ID.d; 00460 I2.d = ID.c + ID.d; 00461 00462 if ( I2.s > 1 ) { 00463 //p2 = p2( 1 / (x+1)); 00464 reverse( I2.p, ID.p); 00465 shift_by_1( I2.p); 00466 } 00467 }else{// midpoint > lowerbound, so perform subdivision like the Descartes method 00468 // Construct the polynomial p(X+midpoint) 00469 // std::cout << "Performing shift by midpoint " << std::endl; 00470 UPOLDSE::shift(I1.p, ID.p, RT(midpoint)); 00471 I1.a = ID.a; 00472 I1.b = ID.a * RT(midpoint) + ID.b; 00473 I1.c = ID.c; 00474 I1.d = ID.c * RT(midpoint) + ID.d; 00475 mult = 0; 00476 // std::cout <<"Transformation is " << I1.a << ":" << I1.b << ":"<< I1.c << ":" << I1.d << std::endl; 00477 00478 if (I1.p[0] == RT(0)) { 00479 RL.push_back( to_FIT( FT( I1.b, I1.d), FT( I1.b, I1.d))); 00480 I1.p /= X; 00481 mult = 1; 00482 } 00483 I1.s = sign_variation( I1.p); 00484 00485 I2.a = ID.b; 00486 I2.b = ID.b + ID.a * RT(midpoint); 00487 I2.c = ID.d; 00488 I2.d = ID.d + ID.c * RT(midpoint); 00489 I2.s = ID.s - I1.s - mult; 00490 00491 // std::cout <<"Transformation is " << I2.a << ":" << I2.b << ":"<< I2.c << ":" << I2.d << std::endl; 00492 // std::cout << "Predicted sign variations of I2 " << ID.s - I1.s - mult << std::endl; 00493 00494 if( I2.s > 1){// Construct the polynomial (X+1)^n p(midpoint/(X+1)) 00495 // std::cout << "Before scaling by " << RT(midpoint) << "poly is "<< ID.p << std::endl; 00496 UPOLDSE::scale( I2.p, ID.p, RT(midpoint)); 00497 // std::cout << "After scaling poly is "<< I2.p << std::endl; 00498 // std::cout << "Sign variations of I2 " << sign_variation(I2.p) << std::endl; 00499 reverse(I2.p); 00500 // std::cout << "poly is " << I2.p << std::endl; 00501 shift_by_1(I2.p); 00502 // std::cout << "After shift poly is " << I2.p << std::endl; 00503 } 00504 }// end of else 00505 00506 if(I2.s > 1){ 00507 if ( I2.p[0] == 0) {//If I2.p[0] = 0 then so is I1.p[0]; so avoid recounting 00508 I2.p /= X; 00509 } 00510 I2.s = sign_variation( I2.p); 00511 } 00512 // std::cout <<"Sign variations " << I1.s << ":" << I2.s << std::endl; 00513 00514 if ( I1.s < I2.s ) { 00515 std::swap( I1, I2); 00516 } 00517 00518 if ( I1.s == 1 ) { 00519 RL.push_back( to_FIT<K>( I1.a, I1.b, I1.c, I1.d, B)); 00520 } else if ( I1.s > 1) { 00521 S.push( IntervalData( I1.a, I1.b, I1.c, I1.d, I1.p, I1.s)); 00522 } 00523 00524 if ( I2.s == 1 ) { 00525 RL.push_back( to_FIT<K>( I2.a, I2.b, I2.c, I2.d, B)); 00526 } else if ( I2.s > 1) { 00527 S.push( IntervalData( I2.a, I2.b, I2.c, I2.d, I2.p, I2.s)); 00528 } 00529 } // while 00530 std::cout << "---------------- iters: " << iters << std::endl; 00531 return; 00532 } 00533 00534 00535 00536 template < typename OutputIterator, typename K > 00537 OutputIterator 00538 CF_solve( const typename K::Poly& f, OutputIterator sol, int mult, K) 00539 { 00540 typedef typename K::Poly Poly; 00541 typedef typename K::FIT FIT; 00542 00543 Seq < FIT > isolating_intervals; 00544 00545 // Check if 0 is a root 00546 int idx; 00547 for (idx = 0; idx <= degree( f); ++idx) { 00548 if ( f[idx] != 0 ) break; 00549 } 00550 00551 Poly p; 00552 if ( idx == 0 ) { p = f; } 00553 else { 00554 p.resize( degree(f) + 1 - idx); 00555 for (int j = idx; j <= degree( f); ++j) 00556 p[j-idx] = f[j]; 00557 } 00558 00559 //std::cout << "Solving: " << p << std::endl; 00560 // Isolate the negative roots 00561 for (int i = 1; i <= degree(p); i+=2) p[i] *= -1; 00562 00563 CF_positive( p, isolating_intervals, false, K()); 00564 00565 // std::cout << "Ok negative" << std::endl; 00566 00567 // Is 0 a root? 00568 if (idx != 0) isolating_intervals.push_back( FIT( 0, 0)); 00569 00570 // Isolate the positive roots 00571 for (int i = 1; i <= degree(p); i+=2) p[i] *= -1; 00572 00573 CF_positive( p, isolating_intervals, true, K()); 00574 00575 // std::cout << "Ok, positive" << std::endl; 00576 00577 // sort( isolating_intervals.begin(), isolating_intervals.end(), CompareFIT()); 00578 00579 // std::cout << "ok sorting" << std::endl; 00580 00581 // for (unsigned i = 0 ; i < isolating_intervals.size(); ++i) { 00582 // if ( singleton( isolating_intervals[i]) ) { 00583 // if ( lower( isolating_intervals[i]) == 0 ) continue; 00584 00585 // Poly t( 2, AsSize()); 00586 // t[1] = denominator( isolating_intervals[i].lower()); 00587 // t[0] = - numerator( isolating_intervals[i].lower()); 00588 00589 // Poly dummy; 00590 // UPOLDSE::div_rem( dummy, p, t); 00591 // p = dummy; 00592 // } 00593 // } 00594 00595 // std::cout << "now p: " << p << ", #nr: " << isolating_intervals.size() << std::endl; 00596 00597 // std::cout << "Done...." << std::endl; 00598 for (unsigned i = 0 ; i < isolating_intervals.size(); ++i) {*sol++ = isolating_intervals[i];} 00599 return sol; 00600 } 00601 00602 00603 00604 00605 00606 template < typename K, 00607 typename OutputIterator > inline 00608 OutputIterator 00609 solve_contfrac( const typename K::Poly& h, OutputIterator sol, K) 00610 { 00611 CF_solve( h, sol, 1, K()); 00612 return sol; 00613 } 00614 00615 00616 00617 } //namespace mmx 00618 00619 00620 #endif // _realroot_SOLVE_CONTFRAC_HPP_ 00621
