solver_uv_continued_fraction.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
 *        P. Dobrowolski, Warsaw University of Technology
 ********************************************************************/
#ifndef realroot_solve_contfrac_hpp
#define realroot_solve_contfrac_hpp

//====================================================================
#include <realroot/GMP.hpp>
#include <realroot/univariate_bounds.hpp>
#include <realroot/sign_variation.hpp>
#include <realroot/Seq.hpp>
#include <realroot/solver.hpp>
#include <realroot/contfrac_intervaldata.hpp>
#include <realroot/contfrac_lowerbound.hpp>

# define TMPL    template<class C, class R, class V>
# define SOLVER  solver<C, ContFrac<R,V> >
//====================================================================
namespace mmx {

template<class M> struct ContFrac{};

//--------------------------------------------------------------------
template<class K, class B, class Extended> struct continued_fraction_isolate;

template<class K, class B> struct continued_fraction_isolate<K, B, false_t> : public K
{
    typedef typename K::integer                                 integer;
    typedef typename K::rational                                rational;
    typedef polynomial<integer, with<MonomialTensor> >          Polynomial;
    typedef typename K::interval_rational                       root_t;
    typedef Seq<root_t>                                         sol_t;
    typedef IntervalData<integer, Polynomial>                   data;

    static inline root_t    as_root(const data& ID)             {return as<root_t>(ID);}
    static inline bool      is_empty(const data& ID)            {return ID.s==0;}
    static inline bool      is_good (const data& ID)            {return ID.s==1;}
    static inline integer   lower_bound(const Polynomial& p)    {return B::lower_bound(p);}
};

template<class K, class B> struct continued_fraction_isolate<K, B, true_t> : public K
{
    typedef typename K::extended_integer                        integer;
    typedef typename K::extended_rational                       rational;
    typedef polynomial<integer, with<MonomialTensor> >          Polynomial;
    typedef typename K::interval_extended_rational              root_t;
    typedef Seq<root_t>                                         sol_t;
    typedef IntervalData<integer, Polynomial>                   data;

    static inline root_t    as_root(const data& ID)             {return as<root_t>(ID);}
    static inline bool      is_empty(const data& ID)            {return ID.s==0;}
    static inline bool      is_good (const data& ID)            {return ID.s==1;}
    static inline integer   lower_bound (const Polynomial& p)   {return B::lower_bound(p);}
};

template<class K, class B, class Extended> struct continued_fraction_approximate;

template<class K, class B> struct continued_fraction_approximate<K, B, false_t> : public K
{
    typedef typename K::integer                                 integer;
    typedef typename K::rational                                rational;
    typedef typename K::interval_rational                       interval_rational;
    typedef polynomial<integer, with<MonomialTensor> >          Polynomial;
    typedef interval_rational                                   root_t;
    typedef Seq<root_t>                                         sol_t;
    typedef IntervalData<integer, Polynomial>                   data;

    static unsigned prec;
    static inline void     set_precision(unsigned n){ prec = n;}
    static inline unsigned get_precision() { return prec; }

    static inline root_t as_root(const data& ID) { return as<root_t>(ID); }

    static inline bool is_empty(const data& ID) {return ID.s==0;}
    static inline bool is_good (const data& ID)
    {
        //      std::cout << __FUNCTION__ <<" "<<ID.p<<" "<<ID.s<<std::endl;
        if(ID.s !=1)
            return false;
        else
        {
            typename data::RT N = ID.a*ID.d-ID.b*ID.c, D=ID.c*ID.d ;
            if(D==0)
                return false;
            if(N<0) N = -N;
            N<<=prec;
            if(D<0) D = -D;
            if(N<D)
                return true;
            else
                return false;
        }
    }

    static inline integer           lower_bound(const Polynomial& p)            {return B::lower_bound(p);}
};

template<class K, class B> struct continued_fraction_approximate<K, B, true_t> : public K
{
    typedef typename K::extended_integer                        integer;
    typedef typename K::extended_rational                       rational;
    typedef typename K::interval_extended_rational              interval_rational;
    typedef polynomial<integer, with<MonomialTensor> >          Polynomial;
    typedef interval_rational                                   root_t;
    typedef Seq<root_t>                                         sol_t;
    typedef IntervalData<integer, Polynomial>                   data;

    static unsigned prec;
    static inline void     set_precision(unsigned n){ prec = n;}
    static inline unsigned get_precision() { return prec; }

    static inline root_t as_root(const data& ID) { return as<root_t>(ID); }

    static inline bool is_empty(const data& ID) {return ID.s==0;}
    static inline bool is_good (const data& ID)
    {
        //      std::cout << __FUNCTION__ <<" "<<ID.p<<" "<<ID.s<<std::endl;
        if(ID.s !=1)
            return false;
        else
        {
            typename data::RT N = ID.a*ID.d-ID.b*ID.c, D=ID.c*ID.d ;
            if(D==0)
                return false;
            if(N<0) N = -N;
            N<<=prec;
            if(D<0) D = -D;
            if(N<D)
                return true;
            else
                return false;
        }
    }

    static inline integer           lower_bound(const Polynomial& p)            {return B::lower_bound(p);}
};

template<class K, class B> unsigned continued_fraction_approximate<K,B,false_t>::prec =
        solver_default_precision;
template<class K, class B> unsigned continued_fraction_approximate<K,B,true_t>::prec =
        solver_default_precision;

template < class K >
struct continued_fraction_subdivision : public K {
    typedef typename K::integer       RT;
    typedef typename K::rational      FT;
    typedef typename K::integer       integer;
    typedef typename K::rational      rational;
    typedef typename K::floating      floating;
    typedef typename K::root_t        root_t;
    typedef typename K::Polynomial    polynomial_integer;
    typedef polynomial< rational, with<MonomialTensor> >  polynomial_rational;
    typedef polynomial< floating, with<MonomialTensor> >  polynomial_floating;
    typedef typename K::data          data;

    static void set_precision(int prec) {K::set_precision(prec);}

    template<class poly>
    static void
    solve_positive(Seq<root_t>& sol, const poly& f, bool posneg) {
        //      std::cout << __FUNCTION__ << " "<<f<<" "<<posneg<<std::endl;

        typedef typename K::data IntervalData;

        IntervalData ID(1, 0, 0, 1, f, 0);
        if (!posneg)
        {
            ID.a=-1;
            for (int i = 1; i <= degree(f); i+=2) ID.p[i] *= -1;
        }
        ID.s = sign_variation(ID.p);

        if ( K::is_empty(ID) ) { return; }
        if ( K::is_good(ID)  ) {
            //std::cout << "A) Sign variation: 1"<<std::endl;;
            FT B = bound_root( f, Cauchy<FT>());
            if ( posneg )
                sol << root_t( FT(0), B);
            else sol << root_t( FT(0), FT(-B));
            return;
        }

        solve_homography(sol,ID, poly());
    }

    template<class output, class data_type, class poly> static void
    solve_homography(output& sol, const data_type& ID, const poly &) {
        //            std::cout << __FUNCTION__ << " "<<ID.p<<std::endl;

        //        typename data_type::Poly X(1,0,1);
        //X[1] = RT(1); X[0]=RT(0);

        int iters  = 0;
        //const RT One(1);
        //FT Zero(0);
        //    std::cout << "Polynomial is " << ID.p << std::endl;

        typedef data_type IntervalData;
        std::stack< IntervalData > S;
        S.push( ID );
        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
            //--------------------------------------------------------------------
            RT a = K::lower_bound(ID.p);
            //    long k = Bd.lower_power_2( ID.p);
            //std::cout<< "\t Lower bound is "<< a<< std::endl;

            if ( a > RT(1) ) {
                scale(ID,a);
                a = RT(1);
            }
            //--------------------------------------------------------------------
            // Step 5 //
            if ( a >= RT(1) ) {
                shift(ID, a);
                //          std::cout<<"Shift by "<<a<<std::endl;

                if ( ID.p[0] == RT(0))
                {
                    div_x(ID.p,1);
                    sol << root_t( to_FT( ID.b, ID.d), to_FT( ID.b, ID.d ));
                }
                ID.s = sign_variation( ID.p);
            }
            //int s= Kernel::init_shift(ID);
            //--------------------------------------------------------------------
            if ( K::is_empty(ID) )
                continue;
            else if( K::is_good(ID) ) {
                sol << K::as_root(ID);
                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;;
                sol << root_t( to_FT( I1.b, I1.d), to_FT(I1.b, I1.d) );
                div_x(I1.p,1);
                r = 1;
            }
            I1.s = sign_variation( I1.p);

            IntervalData I2;
            I2.s = ID.s - I1.s - r;
            reverse_shift_homography(I2,ID);

            // std::cout <<I1.s <<" "<<I2.s<<std::endl;
            // Step 7;
            if ( !K::is_empty(I2) && !K::is_good(I2) ) {

                //p2 = p2( 1 / (x+1));
                reverse( I2.p, ID.p);
                shift_by_1(I2.p);

                if ( I2.p[0] == 0) {
                    div_x(I2.p,1);
                }
                I2.s = sign_variation( I2.p);
            }

            // Step 8
            if ( I1.s < I2.s ) {std::swap( I1, I2); }

            // Step 9
            if ( K::is_good(I1) ) {
                //  std::cout << "C) Sign variation: 1"<<std::endl;;
                sol << K::as_root(I1);
            } else if ( !K::is_empty(I1) ) {
                S.push(I1); //IntervalData( I1.a, I1.b, I1.c, I1.d, I1.p, I1.s));
            }
            // Step 10
            if ( K::is_good(I2) ) {
                //  std::cout << "D) Sign variation: 1"<<std::endl;;
                sol <<  K::as_root(I2);
            } else if ( !K::is_empty(I2) ) {
                S.push(I2); //IntervalData( I2.a, I2.b, I2.c, I2.d, I2.p, I2.s));
            }
        } // while
        //    std::cout << "-------------------- iters: " << iters << std::endl;
        return;
    }

    // solve polynomial_integer or polynomial_extended_integer
    template < class output, class poly> inline static void
    solve_polynomial(output& sol,  const poly& f, int  mult=1) {
        //      std::cout << __FUNCTION__ << " I "<<f<<std::endl;
        Seq < root_t > 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= poly(1,degree(f) - 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;
        solve_positive( isolating_intervals, p, false);
        // Is 0 a root?
        // std::cout << "ok negative" << std::endl;
        if (idx != 0) isolating_intervals << root_t( 0, 0);
        //    for (int i = 1; i <= degree(p); i+=2) p[i] *= -1;
        solve_positive( isolating_intervals, p, true);
        // sort( isolating_intervals.begin(), isolating_intervals.end(), CompareFIT());
        // std::cout << "ok positive" << std::endl;
        // 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 < class output> inline static void
    solve(output& sol,  const polynomial_integer& f) {
        solve_polynomial(sol, f);
    }

    template < class output> inline static void
    solve(output& sol,  const polynomial_rational& f) {
        //      std::cout << __FUNCTION__ << " R "<<f<<std::endl;
        solve_polynomial(sol, univariate::numer<polynomial_integer> (f));
    }

    template < class output> inline static void
    solve(output& sol,  const polynomial_floating& f) {
        //      std::cout << __FUNCTION__ << " R "<<f<<std::endl;
        polynomial_rational p(rational(0),degree(f));

        for(int i=0;i<degree(f)+1;i++) {
            p[i] = as<rational>(f[i]);
        }
        //solve_polynomial(sol,p);
        solve_polynomial(sol, univariate::numer<polynomial_integer> (p));
    }

};

//--------------------------------------------------------------------
template<class C, class Extended> struct solver_approximate_traits;

template<class C>
struct solver_approximate_traits<C, false_t> {
    typedef typename texp::kernelof<C>::T K;
    typedef continued_fraction_subdivision<
    continued_fraction_approximate<K,
    AkritasBound<typename K::integer>, false_t > >  base_t;
    typedef typename K::interval_rational                       interval_t;
};

template<class C>
struct solver_approximate_traits<C, true_t> {
    typedef typename texp::kernelof<C>::T K;
    typedef continued_fraction_subdivision<
    continued_fraction_approximate<K,
    AkritasBound<typename K::extended_integer>, true_t > >  base_t;
    typedef typename K::interval_extended_rational                      interval_t;
};

template<class C>
struct solver<C,ContFrac<Approximate> > {
    typedef typename is_extended<C>::T              Extended;
    typedef solver_approximate_traits<C, Extended>  Traits;
    typedef typename Traits::K                      K;
    typedef typename Traits::base_t                 base_t;
    typedef typename Traits::interval_t             interval_t;
    typedef Seq<interval_t>                         Solutions;
    typedef typename base_t::Polynomial             polynomial_integer;
    typedef typename base_t::polynomial_rational    polynomial_rational;
    typedef typename base_t::polynomial_floating    polynomial_floating;

    static inline void set_precision(unsigned p) {
        continued_fraction_approximate<K,AkritasBound<typename K::integer>,Extended>::set_precision(p);
    }

    template < class output > inline static void
    solve(output& sol, const polynomial_integer& f) {
        base_t::solve(sol, f);
    }

    template < class output > inline
    static void solve(output& sol, const polynomial_rational& f) {
        base_t::solve(sol, f);
    }

    template < class output > inline static void
    solve(output& sol, const polynomial_floating& f) {
        base_t::solve(sol, f);
    }

    template < class output > inline static void
    solve(output& sol, const polynomial_integer& f, int prec) {
        base_t::set_precision(prec);
        solve(sol, f);
        base_t::set_precision(solver_default_precision);
    }

    template <class output> inline static void
    solve(output& sol, const polynomial_integer& f,
          const  typename K::rational& u, const  typename K::rational& v) {
        typedef typename K::integer integer;
        integer
                a= numerator(v),
                b= numerator(u),
                c= denominator(v),
                d= denominator(u);
        IntervalData<integer, polynomial_integer> D =
                as_interval_data(a,b,c,d,f);
        base_t::solve_homography(sol,D,polynomial_integer());
    }

    template <class output> inline static void
    solve(output& sol, const polynomial_rational& f,
          const  typename K::rational& u, const  typename K::rational& v) {
        solve(sol, univariate::numer<polynomial_integer>(f), u, v);
    }


};

//====================================================================
template<class C, class Extended> struct solver_isolate_traits;

template<class C>
struct solver_isolate_traits<C, false_t> {
    typedef typename texp::kernelof<C>::T K;
    typedef continued_fraction_subdivision<
    continued_fraction_isolate<K,
    AkritasBound<typename K::integer>, false_t > >  base_t;
    typedef typename K::interval_rational                       interval_t;
};

template<class C>
struct solver_isolate_traits<C, true_t> {
    typedef typename texp::kernelof<C>::T K;
    typedef continued_fraction_subdivision<
    continued_fraction_isolate<K,
    AkritasBound<typename K::extended_integer>, true_t > >  base_t;
    typedef typename K::interval_extended_rational                      interval_t;
};

template<class C>
struct solver<C,ContFrac<Isolate> > {
    typedef typename is_extended<C>::T          Extended;
    typedef solver_isolate_traits<C,Extended>   Traits;
    typedef typename Traits::base_t             base_t;
    typedef typename Traits::interval_t         interval_t;
    typedef Seq<interval_t>                     Solutions;

    template < class output, class POL> inline static void
    solve(output& sol, const POL& f) {
        base_t::solve(sol,f);
    }
};

//---------------------------------------------------------------------
template<class POL, class M>
typename solver<typename POL::Scalar, ContFrac<M> >::Solutions
solve (const POL& p, const ContFrac<M>& slv) {
    typedef typename POL::Scalar                             Scalar;
    typedef solver<Scalar, ContFrac<M> >                     Solver;
    typedef typename solver<Scalar, ContFrac<M> >::Solutions Solutions;
    Solutions sol;
    Solver::solve(sol,p);
    return sol;
}

//---------------------------------------------------------------------
template<class POL, class M, class B>
typename solver<typename POL::Scalar, ContFrac<M> >::Solutions
solve (const POL& p, const ContFrac<M>& slv, const B& b1, const B& b2) {
    typedef solver<B, ContFrac<M> >                     Solver;
    typedef typename solver<B, ContFrac<M> >::Solutions Solutions;
    Solutions sol;
    Solver::solve(sol,p,b1,b2);
    return sol;
}
//--------------------------------------------------------------------

} //namespace mmx
//====================================================================
# undef TMPL
# undef SOLVER
//====================================================================
#endif // _realroot_solve_contfrac_hpp_