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

Go to the documentation of this file.

#ifndef synaps_usolve_contfrac_IntervalData_h
#define synaps_usolve_contfrac_IntervalData_h

#include <realroot/Interval.hpp>
#include <realroot/univariate.hpp>
namespace mmx {

    template<class POL> void div_x(POL& p, int k)
    {
      for(int i=0;i< int(p.size()-k);i++)
        p[i]=p[i+k];
      for(int i=p.size()-k; i<(int)p.size();i++)
        p[i]=0; 
    }

    template < typename RT, typename Poly >
    struct IntervalData 
    {
        RT a, b, c, d;
        Poly  p;
        unsigned long s;
  
        IntervalData( const RT& a_, const RT& b_, const RT& c_, const RT& d_,
                      const Poly& p_, unsigned long s_)
            : a(a_), b(b_), c(c_), d(d_), p(p_), s(s_) {}
  
        IntervalData() {}
  
    };
  
    template < typename FT > inline
    Interval< FT > as_interval( const FT& a, const FT& b)
    {
        if ( a <= b ) 
            return Interval<FT>( a, b);
        return Interval<FT>( b, a);
    }

    template < typename RT, typename Poly >
    IntervalData<RT,Poly> as_interval_data(const RT& a, const RT& b, const RT&c, const RT&d,
                                           const Poly& f)
    {

      Poly  X = Poly("x")*a+Poly(b), Y=Poly("x")*c+Poly(d);

      Poly  F = univariate::eval_homogeneous(f,X,Y);
      return IntervalData<RT,Poly> (a,b,c,d,F,degree(f));
    } 

#ifndef WITH_AS
#define WITH_AS
    template<typename T,typename F>
    struct as_helper {
        static inline T cv (const F& x) { return x; }
    };

    template<typename T,typename F> inline T
    as (const F& x) {  return as_helper<T,F>::cv (x); }
#endif

  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); 
    }
  };






  // template <> inline rational as<rational>(const integer& a, const integer& b)
  //  {
  //    return rational(a)/b;
  //  }

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



  template < typename FT >
  Interval< FT > inline
  to_FIT( const FT& a, const FT& b)
  {
    if ( a <= b ) 
      return Interval<FT>( a, b);
    return Interval<FT>( b, a);
  }


  template < typename K >
  typename K::FIT
  to_FIT( const typename K::RT& a, const typename K::RT& b, 
          const typename K::RT& c, const typename K::RT& d, 
          const typename K::FT& B)
  {
    typedef typename K::RT    RT;
    typedef typename K::FT    FT;
    typedef typename K::FIT   FIT;
    
    
    FT left;
    FT right;
    
    if ( c == RT(0)) left = sign(a) * B;
    //  else left = as<FT>(a, c);
    else left = to_FT(a,c);
  
    if ( d == RT(0)) right = sign(b) * B;
        //  else right = as<FT>(b, d);
    else right = to_FT(b,d);
  
    return to_FIT( left, right);
  }


  template < typename Poly> inline 
  void 
  reverse(Poly & R, const Poly& P) { 
    R=P;  reverse(R);
  }
  
  /* Replace R by its reverse */
  template < typename Poly> inline 
  void 
  reverse(Poly & R)
  { 
    int deg = degree(R);
    Poly temp(R);
    for (int i = 0; i <= deg; ++i)
      R[i] = temp[deg-i];
  }
  
  /* Replaces P by the polynomial P(X+1) */
  template < typename RT, typename Poly > inline
  void reverse( IntervalData<RT,Poly>& I1)
  {
    reverse(I1.p);
    std::swap(I1.a,I1.b);
    std::swap(I1.c,I1.d);
    return; 
  }



  /* Replaces P by the polynomial P(X+1) */
  template < typename Poly > inline
  void shift_by_1( Poly& R, const Poly& P)
  {
    R=P;
    int deg = degree(P);
    for (int j = deg-1; j >= 0; j--) {
      for (int i = j; i < deg; i++) {
        R[i] +=  R[i+1];
      }

    }
    
    return;
  }
  
  /* Replaces P by the polynomial P(X+1) */
  template < typename RT, typename Poly > inline
  void shift_by_1( IntervalData<RT,Poly>& I1, const IntervalData<RT,Poly>& I2)
  {
    shift_by_1(I1.p,I2.p);
    I1.a = I2.a;
    I1.b = I2.a + I2.b;
    I1.c = I2.c;
    I1.d = I2.c + I2.d;
    I1.s = sign_variation(I1.p);
    return; 
  }

  /* Replaces P by the polynomial P(X+1) */
  template < typename RT, typename Poly > inline
  void shift_by_1( IntervalData<RT,Poly>& I1)
  {
    shift_by_1(I1.p);
    I1.b += I1.a;
    I1.d += I1.c;
    I1.s = sign_variation(I1.p);
  }
  
  /* Replaces P by the polynomial P(X+1) */
  template < typename Poly > inline
  void shift_by_1(Poly& R)
    
  {
    long deg = degree( R);
    for (int j = deg-1; j >= 0; j--) {
      for (int i = j; i < deg; i++) {
        R[i] +=  R[i+1]; 
      }
    }
    return; 
  }
  

  template < typename RT, typename Poly > inline
  void scale( IntervalData<RT,Poly>& ID, const RT& a)
  {
    univariate::scale( ID.p, ID.p, a);
    ID.a *= a;
    ID.c *= a;
  }
  
  
  template < typename RT, typename Poly > inline
  void scale( IntervalData<RT,Poly>& ID, long k)
  {
    RT a = pow( RT( 2), k);
    
    univariate::scale( ID.p, ID.p, a);
    ID.a *= a;
    ID.c *= a;
  }
  
  
  template < typename RT, typename Poly > inline
  void inv_scale( IntervalData<RT,Poly>& ID, const RT& a)
  {
    univariate::inv_scale( ID.p, ID.p, a);
    ID.b *= a;
    ID.d *= a;
  }
  
  
  template < typename RT, typename Poly > inline
  void shift( IntervalData<RT,Poly>& ID, const RT& a)
  {
    univariate::shift( ID.p, ID.p, a);
    ID.s = sign_variation(ID.p);
    ID.b += a * ID.a;
    ID.d += ID.c * a;
  }
  
  
  /* Replaces P by the polynomial (X+1)^d*P(1/(X+1)) */
  template < typename RT, typename Poly > inline
  void reverse_shift_homography( IntervalData<RT,Poly>& I2, const IntervalData<RT,Poly>& ID)
  {
    I2.a = ID.b;
    I2.b = ID.a + ID.b;
    I2.c = ID.d;
    I2.d = ID.c + ID.d;
  }
  
  /* Replaces P by the polynomial (X+1)^d*P(1/(X+1)) */
  template < typename RT, typename Poly > inline
  void reverse_shift( IntervalData<RT,Poly>& I2, const IntervalData<RT,Poly>& ID)
  {
    I2.a = ID.b;
    I2.b = ID.a + ID.b;
    I2.c = ID.d;
    I2.d = ID.c + ID.d;
    I2.p= ID.p;
    reverse(I2);
    shift_by_1(I2);
    I2.s = sign_variation(I2.p);
  }

} //namespace mmx
#endif //synaps_usolve_contfrac_IntervalData_h