realroot: include/realroot/contfrac.hpp Source File

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