contfrac.hpp

Go to the documentation of this file.

/*******************************************************************
 *   This file is part of the source code of the realroot kernel.
 *   Author(s): 
 *        E. Tsigaridas <et@di.uoa.gr>
 *           Department of Informatics and Telecommunications
 *           University of Athens (Greece).
 *   Modification:
 *        B. Mourrain, GALAAD, INRIA
 *        V. Sharma, GALAAD, INRIA
 ********************************************************************/


#ifndef _realroot_SOLVE_CONTFRAC_HPP_
#define _realroot_SOLVE_CONTFRAC_HPP_

#include <realroot/Interval.hpp>
//#include <realroot/traits.hpp>
#include <realroot/upoldse.hpp>
#include <realroot/upoldse_bound.hpp>
#include <realroot/sign.hpp>
#include <realroot/Seq.hpp>

#include <realroot/contfrac_intervaldata.hpp>
#include <realroot/contfrac_lowerbound.hpp>



namespace mmx {


    namespace meta {
        
        class Kioustelidis_bound_1
        {

          public:
            Kioustelidis_bound_1() {}


            template <typename NT>
            long upper( const upoldse<NT>& p)
            {
                return foo( p.begin(), p.end());
                // return foo( begin( p), end( p));
            }

            template <typename NT>
            long lower( const upoldse<NT>& p)
            {
                return -foo( p.rbegin(), p.rend());
            }
    
            template <class IT>
            long foo( IT first, IT last)
            {
                typedef typename std::iterator_traits<IT>::difference_type     size_type;
        
                size_type deg = std::distance( first, last);
                assert( deg > 0 );

                int s = sign( *(last-1));
                long lan = bitsize( *(last-1));

                bool first_time = true;
                long maxpow = 0;

                int i = 1;
                for ( IT it = (last-2); it != first-1; --it, ++i) {
                    if ( s * sign( *it) < 0) {
                        long p = bitsize( abs( *it)) - lan - 1;
                        long q = p / i;
                        if ( first_time ) {
                            first_time = false;
                            maxpow = q + 2;
                        } else {
                            maxpow = std::max( maxpow, q + 2);
                        }
                    }
                }
                return maxpow;
            }
        };

        template <typename T> 
        void times_2_to_k( T& r, const T& a, long k = 1)
        {
            r = a << k;
            
        }

        template <typename T>
        unsigned long int bitsize( const T& a)
        {
            return std::log( static_cast<double>(a));
        }        

#define __GMP_EXPR ::__gmp_expr<T, U>
#define TMPL template <class T, class U> inline
        
        TMPL 
        mpz_class times_2_to_k( const __GMP_EXPR& a, long k = 1)
        {
            mpz_class b;
            mpz_mul_2exp( b.get_mpz_t(), a.get_mpz_t(), k );
            return b;
        }
        
        TMPL 
        void times_2_to_k( __GMP_EXPR& r,  const __GMP_EXPR& a, long k = 1)
        {
            mpz_mul_2exp( r.get_mpz_t(), a.get_mpz_t(), k );
        }

        TMPL
        unsigned long int bitsize( const __GMP_EXPR& a) 
        {
            return mpz_sizeinbase( mpz_class( a).get_mpz_t(), 2 );
        }

#undef TMPL
#undef __GMP_EXPR

        
        
        template < typename RT, typename Poly > inline
        void do_scale( IntervalData<RT,Poly>& ID, long k)
        {
            assert( k > 0 );
            
            for (int i = 0; i <= degree( ID.p); ++i) {
                times_2_to_k( ID.p[i], ID.p[i], k*i);
            }
            times_2_to_k( ID.a, ID.a, k);
            times_2_to_k( ID.c, ID.c, k);
        }

    } //meta
    

    template<class C> typename meta::rationalof<C>::T to_rational(const C& a, const C& b)
    {
        typedef typename meta::rationalof<C>::T rational;
        return  (rational(a)/b);
    }

    template < class NT , class LowerBound = AkritasBound<NT> >
    struct ContFrac_t 
    {
        typedef NT                 RT;
        typedef typename meta::rationalof<RT>::T     FT;
        typedef Interval<FT>       FIT;
        typedef Interval<FT>       root_t;
        typedef upoldse<NT>        Poly;
        typedef LowerBound         bound;
        typedef ContFrac_t         self_t;
    };


    template < class NT, class LB > inline
    Seq< typename ContFrac<NT,LB>::root_t > 
    solve( const typename ContFrac<NT>::Poly& f, ContFrac<NT,LB>)
    {
        typedef ContFrac<NT,LB>        K;
        typedef typename K::root_t     root_t;
    
        Seq<root_t> sol;
        //    std::cout << "**In this" << std::endl; 
        solve_contfrac(f, std::back_inserter(sol.rep()), K());
        //  std::stable_sort( sol.begin(), sol.end(), ALGEBRAIC::Refine_compare());
    
        return sol;
    }


    template < typename K >
    void
    CF_positive( const typename K::Poly& f, Seq< typename K::FIT>& RL, bool posneg, K)
    {
        //  std::cout << __FUNCTION__ << std::endl; 
        typedef typename K::RT    RT;
        typedef typename K::FT    FT;
        typedef typename K::FT    rational_t;
        typedef typename K::FIT   FIT;
        typedef typename K::Poly  Poly;
  
        typedef IntervalData< RT, Poly> IntervalData;
        
        meta ::Kioustelidis_bound_1 bound;
        
        
        int iters  = 0;
        const RT One( 1);
        FT Zero(0);
 
        // FT B = bound_root( f, Cauchy<RT>());
        
        FT B = meta::times_2_to_k( RT( 1), bound.upper( f));
        
        
        unsigned long s = sign_variation(f);

        if ( s == 0 ) { return; }
        if ( s == 1 ) {
            //      std::cout << "A) Sign variation: 1"<<std::endl;;
            if ( posneg )
            {
                //    std::cout<<to_FIT( FT(0), B)<<std::endl;;
                RL.push_back( to_FIT( FT(0), B));
            }
            else RL.push_back( to_FIT( FT(0), FT(-B)));
            return;
        }
        Poly X(2, AsSize()); X[1] = RT(1); X[0]=RT(0);

        std::stack< IntervalData > S;
        if ( posneg) S.push( IntervalData( 1, 0, 0, 1, f, s));
        else S.push( IntervalData( -1, 0, 0, 1, f, s));

        while ( !S.empty() ) {
            ++iters;
            IntervalData ID = S.top();
            S.pop();

            //std::cout << "Polynomial is " << ID.p << std::endl;
            // Steps 3 - 4: Choose the method for computing the lower bound 


            //typename K::bound Bd;
            // Cauchy<RT> Bd;
            //HongBound<RT> Bd;
      

            //    RT M = Bd.upper(ID.p);
            //    scale(ID,M);
            //    reverse(ID);
            //    translate_by_1(ID);
            //    reverse(ID);

            //     Interval<rational_t> q=Bd.range(ID.p);
            //     std::cout<<"Sign var: "<< sign_variation(ID.p)<<" ";
            //     RT a, b;
            //     if(q<0)
            //       {
            //      continue;
            //       }
            //     else if(q!=0)
            //       {
            //      a=numerator(q.lower());
            //      b=denominator(q.lower());
            //      inv_scale(ID,b);
            //      shift(ID, a);
            //      scale(ID,b);
            //  let::assign(a,rational_t(q.upper()-q.lower()));
            //  shift(ID,a);
            //      a=One;
            //       }

            //    std::cout<<"("<< sign_variation(ID.p)<<")"<<std::endl;;


            // RT a = Bd.lower(ID.p);
            //long k = bound.lower( ID.p);
            long k = AkritasBound<RT>().lower_power_2( ID.p);
            // std::cout << "k: " << k << std::endl; 
            //    std::cout<< "\t Lower bound is "<< a<< std::endl;
            
            if (k > 0) {
                meta::do_scale( ID, k);
            }
            
            // if ( a > One ) { scale(ID,a); a = One; }
    
            // Step 5 //
            // if ( a >= One ) 
            if ( k >= 0 ) {
                shift_by_1(ID);
                //shift(ID, a);
      
                if ( ID.p[0] == RT(0)) {
                    RL.push_back( to_FIT( to_rational( ID.b, ID.d), to_rational( ID.b, ID.d )));
                    ID.p /= X;
                }
                ID.s = sign_variation( ID.p);
                if ( ID.s == 0 ) { continue; }
                if ( ID.s == 1 ) {
                    //    std::cout << "B) Sign variation: 1"<<std::endl;
                    RL.push_back( to_FIT<K>( ID.a, ID.b, ID.c, ID.d, B));
                    continue;
                }
            }
    
            // Step 6
            IntervalData I1;
            shift_by_1( I1, ID);

            unsigned long r  = 0;
    
            if (I1.p[0] == RT(0)) {
                //  std::cout << "Zero at end point"<<std::endl;;
                RL.push_back( to_FIT( to_rational( I1.b, I1.d), to_rational(I1.b, I1.d) ));
                I1.p /= X;
                r = 1;
            }
            I1.s = sign_variation( I1.p);

            IntervalData I2;
        
            I2.s = ID.s - I1.s - r;

            I2.a = ID.b;
            I2.b = ID.a + ID.b;
            I2.c = ID.d;
            I2.d = ID.c + ID.d;
    
            // Step 7;
            if ( I2.s > 1 ) {
                //p2 = p2( 1 / (x+1));
                reverse( I2.p, ID.p);
                shift_by_1( I2.p);

                if ( I2.p[0] == 0) {
                    I2.p /= X;
                }
                I2.s = sign_variation( I2.p);
            }

            // Step 8
            if ( I1.s < I2.s ) {
                std::swap( I1, I2);
            }
    
            // Step 9
            if ( I1.s == 1 ) {
                //  std::cout << "C) Sign variation: 1"<<std::endl;;
                RL.push_back( to_FIT<K>( I1.a, I1.b, I1.c, I1.d, B));
            } else if ( I1.s > 1) {
                S.push( IntervalData( I1.a, I1.b, I1.c, I1.d, I1.p, I1.s));
            }

            // Step 10
            if ( I2.s == 1 ) {
                //  std::cout << "D) Sign variation: 1"<<std::endl;;
                RL.push_back( to_FIT<K>( I2.a, I2.b, I2.c, I2.d, B));
            } else if ( I2.s > 1) {
                S.push( IntervalData( I2.a, I2.b, I2.c, I2.d, I2.p, I2.s));
            }
        } // while
        //    std::cout << "-------------------- iters: " << iters << std::endl; 
        return;  
    }

    // A continued fraction algorithm that combines the Descartes method and
    // Akritas' algorithm. The basic idea is to take the better of performing
    // a Taylor shift or subdividing into half.
    template < typename K >
    void
    MCF_positive( const typename K::Poly& f, Seq< typename K::root_t>& RL, bool posneg, K)
    {
        //  std::cout << __FUNCTION__ << std::endl; 
        typedef typename K::RT    RT;
        typedef typename K::FT    FT;
        typedef typename K::root_t   FIT;
        typedef typename K::Poly  Poly;
  
        typedef IntervalData< RT, Poly> IntervalData;
  
        int iters  = 0;
        //std::cout <<"Polynomial is "<< f << std::endl;

        FT B = bound_root( f, Cauchy<FT>());
        //  std::cout << "Upper bound is "<< B << " Positive or negative " << posneg << std::endl;
        unsigned long s = sign_variation( f);
        if ( s == 0 ) { return; }
        if ( s == 1 ) {
            if ( posneg ) RL.push_back( to_FIT( FT(0), B));
            else RL.push_back( to_FIT( FT(0), FT(-B)));
            return;
        }

        Poly X(2, AsSize()); X[1] = RT(1);

        std::stack< IntervalData > S;
        if ( posneg) S.push( IntervalData( 1, 0, 0, 1, f, s));
        else S.push( IntervalData( -1, 0, 0, 1, f, s));

        IntervalData ID, I1, I2;// ID is the data associated with the parent and
        // I1, I2 are the data associated with its children

        RT lowerBound;
        FT midpoint; // The point corresponding to the midpoint of the interval 
        // with endpoints (a/c), (b/d).
        FT temp; // A temporary variable 
        unsigned long mult = 0; // multiplicity of an exact root.

        while ( !S.empty() ) {
            ++iters;
            ID = S.top();
            S.pop();
            //    std::cout <<"Transformation is " << ID.a << ":" << ID.b << ":"<< ID.c << ":" << ID.d << std::endl;
    
            // Steps 3 - 4
            lowerBound =  CauchyMinimumBound( ID.p);
            std::cout << "lower bound is "<< lowerBound << std::endl;
    
            if(ID.c == 0){// then ID.a is not equal to zero
                midpoint = FT(ID.d * B - ID.b*sign(ID.a), 2* abs(ID.a));
                //std::cout << "Midpoint = " << midpoint << std::endl;
            }else if (ID.d > ID.c ){
                midpoint = FT(ID.d, ID.c);
            }else{}

            if(midpoint < 1)  midpoint = 1;
            std::cout << "Midpoint = " << midpoint << std::endl;
    
            if(lowerBound >= midpoint){
                if ( lowerBound > 1 ) {
                    UPOLDSE::scale( ID.p, ID.p, lowerBound);
                    ID.a *= lowerBound;
                    ID.c *= lowerBound;
                    lowerBound = 1;
                }
    
                if ( lowerBound == 1 ) {
                    shift_by_1( ID.p, ID.p);// If we scaled above then this is required
                    ID.b += ID.a;               // to bring the origin to lowerBound
                    ID.d += ID.c;
      
                    if ( ID.p[0] == RT(0)) {
                        RL.push_back( to_FIT( FT( ID.b, ID.d), FT( ID.b, ID.d)));
                        ID.p /= X;
                    }
                    ID.s = sign_variation( ID.p);
                    if ( ID.s == 0 ) { continue; }
                    if ( ID.s == 1 ) {
                        RL.push_back( to_FIT<K>( ID.a, ID.b, ID.c, ID.d, B));
                        continue;
                    }
                }

                shift_by_1( I1.p, ID.p);
                I1.a = ID.a;
                I1.b = ID.a + ID.b;
                I1.c = ID.c;
                I1.d = ID.c + ID.d;
    
                if (I1.p[0] == RT(0)) {
                    RL.push_back( to_FIT( FT( I1.b, I1.d), FT( I1.b, I1.d)));
                    I1.p /= X;
                    mult = 1;
                }
                I1.s = sign_variation( I1.p);

                I2.s = ID.s - I1.s - mult;
                I2.a = ID.b;
                I2.b = ID.a + ID.b;
                I2.c = ID.d;
                I2.d = ID.c + ID.d;
    
                if ( I2.s > 1 ) {
                    //p2 = p2( 1 / (x+1));
                    reverse( I2.p, ID.p);
                    shift_by_1( I2.p);
                }
            }else{// midpoint > lowerbound, so perform subdivision like the Descartes method
                // Construct the polynomial p(X+midpoint)
                //      std::cout << "Performing shift by midpoint " << std::endl;
                UPOLDSE::shift(I1.p, ID.p, RT(midpoint));
                I1.a = ID.a;
                I1.b = ID.a * RT(midpoint) + ID.b;
                I1.c = ID.c;
                I1.d = ID.c * RT(midpoint) + ID.d;
                mult = 0;
                //      std::cout <<"Transformation is " << I1.a << ":" << I1.b << ":"<< I1.c << ":" << I1.d << std::endl;

                if (I1.p[0] == RT(0)) {
                    RL.push_back( to_FIT( FT( I1.b, I1.d), FT( I1.b, I1.d)));
                    I1.p /= X;
                    mult = 1;
                }
                I1.s = sign_variation( I1.p);

                I2.a = ID.b;
                I2.b = ID.b + ID.a * RT(midpoint);
                I2.c = ID.d;
                I2.d = ID.d + ID.c * RT(midpoint);
                I2.s = ID.s - I1.s - mult;

                //      std::cout <<"Transformation is " << I2.a << ":" << I2.b << ":"<< I2.c << ":" << I2.d << std::endl;
                //      std::cout << "Predicted sign variations of I2 " << ID.s - I1.s - mult << std::endl;
      
                if( I2.s > 1){// Construct the polynomial (X+1)^n p(midpoint/(X+1))
                    //  std::cout << "Before scaling by " << RT(midpoint) << "poly is "<< ID.p << std::endl;
                    UPOLDSE::scale( I2.p, ID.p, RT(midpoint));
                    //  std::cout << "After scaling poly is "<< I2.p << std::endl;  
                    //  std::cout << "Sign variations of I2 " << sign_variation(I2.p) << std::endl;
                    reverse(I2.p);
                    //  std::cout << "poly is " << I2.p << std::endl;
                    shift_by_1(I2.p);
                    //  std::cout << "After shift poly is " << I2.p << std::endl;
                }
            }// end of else

            if(I2.s > 1){
                if ( I2.p[0] == 0) {//If I2.p[0] = 0 then so is I1.p[0]; so avoid recounting
                    I2.p /= X;
                }
                I2.s = sign_variation( I2.p);
            }
            //    std::cout <<"Sign variations " << I1.s << ":" << I2.s << std::endl;

            if ( I1.s < I2.s ) {
                std::swap( I1, I2);
            }
    
            if ( I1.s == 1 ) {
                RL.push_back( to_FIT<K>( I1.a, I1.b, I1.c, I1.d, B));
            } else if ( I1.s > 1) {
                S.push( IntervalData( I1.a, I1.b, I1.c, I1.d, I1.p, I1.s));
            }
    
            if ( I2.s == 1 ) {
                RL.push_back( to_FIT<K>( I2.a, I2.b, I2.c, I2.d, B));
            } else if ( I2.s > 1) {
                S.push( IntervalData( I2.a, I2.b, I2.c, I2.d, I2.p, I2.s));
            }
        } // while
        std::cout << "---------------- iters: " << iters << std::endl; 
        return;  
    }



    template < typename OutputIterator, typename K >
    OutputIterator
    CF_solve( const typename K::Poly& f, OutputIterator sol, int  mult, K)
    {
        typedef typename K::Poly  Poly;
        typedef typename K::FIT FIT;

        Seq < FIT > isolating_intervals;
  
        // Check if 0 is a root
        int idx;
        for (idx = 0; idx <= degree( f); ++idx) {
            if ( f[idx] != 0 ) break;
        }
  
        Poly p;
        if ( idx == 0 ) { p = f; } 
        else {
            p.resize( degree(f) + 1 - idx);
            for (int j = idx; j <= degree( f); ++j) 
                p[j-idx] = f[j];
        }
  
        //std::cout << "Solving: " << p << std::endl; 
        // Isolate the negative roots
        for (int i = 1; i <= degree(p); i+=2) p[i] *= -1;

        CF_positive( p, isolating_intervals, false, K());

        //  std::cout << "Ok negative" << std::endl; 

        // Is 0 a root?
        if (idx != 0) isolating_intervals.push_back( FIT( 0, 0));
  
        // Isolate the positive roots
        for (int i = 1; i <= degree(p); i+=2) p[i] *= -1;

        CF_positive( p, isolating_intervals, true, K());

        //  std::cout << "Ok, positive" << std::endl; 

        //  sort( isolating_intervals.begin(), isolating_intervals.end(), CompareFIT());

        //  std::cout << "ok sorting" << std::endl; 

  //       for (unsigned i = 0 ; i < isolating_intervals.size(); ++i) {
//             if ( singleton( isolating_intervals[i]) ) {
//                 if ( lower( isolating_intervals[i]) == 0 ) continue;
      
//                 Poly t( 2, AsSize());
//                 t[1] = denominator( isolating_intervals[i].lower());
//                 t[0] = - numerator( isolating_intervals[i].lower());
      
//                 Poly dummy;
//                 UPOLDSE::div_rem( dummy, p, t);
//                 p = dummy;
//             }
//         }
  
        //  std::cout << "now p: " << p << ",  #nr: " << isolating_intervals.size() << std::endl; 

        //  std::cout << "Done...." << std::endl; 
        for (unsigned i = 0 ; i < isolating_intervals.size(); ++i) {*sol++ = isolating_intervals[i];}
        return sol;
    }





    template < typename K,
               typename OutputIterator > inline
    OutputIterator
    solve_contfrac( const typename K::Poly& h, OutputIterator sol, K)
    {
        CF_solve( h, sol, 1, K()); 
        return sol;
    }



} //namespace mmx 


#endif // _realroot_SOLVE_CONTFRAC_HPP_