/Users/mourrain/Devel/mmx/realroot/include/realroot/solver_univariate_glue.hpp

Go to the documentation of this file.

#include <basix/glue.hpp>
#include <basix/vector.hpp>
#include <basix/vector_sort.hpp>
#include <numerix/kernel.hpp>
#include <realroot/Seq.hpp>
#include <realroot/Interval_glue.hpp>
#include <realroot/polynomial_tensor.hpp>
#include <realroot/solver_uv_sleeve.hpp>
#include <realroot/solver_uv_continued_fraction.hpp>
#include <realroot/solver_uv_interval_newton.hpp>


#define TMPL        template<class C> 
#define RING        ring<C,MonomialTensor>
#define POLYNOMIAL  polynomial< C, with<MonomialTensor> >

namespace mmx {

  template<typename FT, typename RT, typename Poly> 
  struct as_helper<interval<FT>,IntervalData<RT,Poly> > {
    static inline interval<FT>
    cv (const IntervalData<RT,Poly>& I) { 
      FT left, right;
      if ( I.c == RT(0)) 
        left = sign(I.a) * bound_root( I.p, Cauchy<FT>());
      else 
        left = FT(I.a)/I.c;
      
      if ( I.d == RT(0)) 
        right = sign(I.b) * bound_root( I.p, Cauchy<FT>()); 
      else 
        right = FT(I.b)/I.d;
      //      std::cout << I.a<<"/"<<I.c<<", "<<I.b<<"/"<<I.d<<std::endl;
      return interval<FT>(left, right); 
    }
  };

  TMPL vector<generic> 
  solver_univariate_contfrac_prec(const POLYNOMIAL & p, const int& prec) {
    typedef Numerix                             kernel;
    typedef typename kernel::integer            Integer;
    typedef typename kernel::rational           Rational;
    typedef typename kernel::interval_rational  interval_rational;
    typedef solver< Integer, ContFrac<Approximate> > Solver;

    vector<interval_rational> sol;
    Solver::set_precision(prec);
    Solver::solve(sol,p);
    sort(sol);

    vector<generic> r;
    for(unsigned i=0;i<sol.size();i++) 
      r <<as<generic>(interval< Rational > (lower(sol[i]),upper(sol[i])));
    //r <<as<generic>(sol[i]);
    return r;
  }

  //--------------------------------------------------------------------
  TMPL vector<generic> solver_univariate_contfrac(const POLYNOMIAL & p) {
    return solver_univariate_contfrac_prec(p,mmx_bit_precision);
  }
  
  //--------------------------------------------------------------------
  TMPL vector<generic> 
  solver_univariate_contfrac(const POLYNOMIAL & p, const interval<rational>& D) {
    typedef Numerix                        kernel;
    typedef typename kernel::integer       Integer;
    typedef typename kernel::rational      Rational;
    typedef typename kernel::interval_rational interval_rational;
    typedef solver<Integer, ContFrac<Approximate> >       Solver;

    vector<interval_rational> sol;
    Solver::set_precision(mmx_bit_precision);
    Solver::solve(sol,p,lower(D),upper(D));
    sort(sol);
    vector<generic> r;
    for(unsigned i=0;i<N(sol);i++) 
      r <<as<generic>(interval< rational > (lower(sol[i]),upper(sol[i])));
    return r;
  }
  //--------------------------------------------------------------------


  TMPL vector<generic> 
  solver_univariate_contfrac_approx_prec(const POLYNOMIAL & p, const int& prec) {
    typedef Numerix                            kernel;
    typedef typename kernel::integer           Integer;
    typedef typename kernel::rational          rational;
    typedef typename kernel::interval_rational interval_rational;
    typedef solver<Integer,ContFrac<Approximate> >  Solver;


    Solver::set_precision(prec);

    vector<interval_rational> sol;
    Solver::solve(sol,p);
    sort(sol);
    long int old_prec = mmx_bit_precision;
    vector<generic> r;
    mmx_bit_precision = prec;
    for(unsigned i=0;i<sol.size();i++) {
      r << as<generic>(as<floating<> >((lower(sol[i])+upper(sol[i]))/2));
    }
    mmx_bit_precision = old_prec;
    return r;
  }
  //--------------------------------------------------------------------
  TMPL vector<generic> 
  solver_univariate_contfrac_approx(const POLYNOMIAL & p) {
    return solver_univariate_contfrac_approx_prec(p,mmx_bit_precision);
  }
  //--------------------------------------------------------------------
  TMPL vector<generic> 
  solver_univariate_sleeve(const POLYNOMIAL & p) {
    typedef Numerix                   kernel;
    typedef typename kernel::integer  Integer;
    typedef typename kernel::rational rational;
    typedef typename kernel::interval_rational interval_rational;
    typedef kernel::interval_rational interval_rational;
    Seq<interval_rational> sol;
    solver<RING,Sleeve<Isolate> >::solve(sol,p);
    vector<generic> r;
    for(unsigned i=0;i<sol.size();i++) r << as<generic>(sol[i]);
    return r;
  }
  //--------------------------------------------------------------------
  TMPL interval<floating<> > 
  solver_univariate_newton(const POLYNOMIAL & p, const interval<floating<> >& I) {
    typedef interval<floating<> > interval_t;
    IntervalNewton<interval_t,C> Nwt(I);
    Seq<interval_t> sol =solve(p,Nwt);
    
    return sol[0];
  }

} // namespace mmx
//--------------------------------------------------------------------
#undef TMPL
#undef RING
#undef POLYNOMIAL