Interval.hpp

Go to the documentation of this file.

#ifndef mmx_interval_hpp
#define mmx_interval_hpp

#include <cassert>
#include <realroot/print.hpp>
#include <realroot/texp.hpp>
#include <realroot/texp_expression.hpp>
#include <realroot/rounding_mode.hpp>

namespace mmx {

  using namespace texp;

namespace  numerics { template<class T, int r> struct interval_base; }

template<class T, int r = 3 >
struct Interval : numerics::interval_base<T,((unsigned)r)%4> {
    typedef T value_type;
    typedef T coeff_t;
    typedef T boundary_type;
    typedef numerics::interval_base<T,((unsigned)r)%4> base_t;
    typedef typename base_t::rnd rnd_t;
    typedef Interval<T,r> self_t;
    struct extended { T a,b; extended( const T& c, const T& d ) { a = c; b = d; }; };
    T m, M;
    template<class S> Interval( const texp::template_expression<S>& e );
    Interval();    
    Interval( int n );
    Interval( unsigned n );
    Interval( const T& x );
    Interval( const T& a, const T& b );
    Interval( const char * s );
    //Interval& operator=( const char * s ){ let::assign(*this,s); return *this; }
    Interval& operator=( const Interval & i )
    {
      m = i.m;
      M = i.M;
      return *this;
    };
    
    template<class X>
    Interval & operator=( const texp::template_expression<X> & x ) ;
    
    void define( const T& m, const T& M ) { this->m = m; this->M = M; };
    void assign( const T& m, const T& M ) { define(m,M); };
    void assign( const Interval& b ) { *this = b; };

    template<class X> bool operator==( const X& k) const { return m == M && M == k ; };
    template<class X, int r_> bool operator==( const Interval<X,r_> & b ) const 
    { return b.m == m && b.M == M ; };
    template<class X> inline bool operator>( const X& x ) const { return m > x; };
    template<class X> inline bool operator>=( const X& x ) const { return  m >= x; };
    template<class X> inline bool operator<( const X& x ) const { return M < x; };
    template<class X> inline bool operator<=( const X& x ) const { return M <= x; };
    template<class X> inline bool operator!=(const X& x ) const { return !(*this == x); };
    template<class X, int _r> bool operator<( const Interval<X,_r>& i ) const { return M < i.m; };
    template<class X, int _r> bool operator>( const Interval<X,_r>& i ) const { return m > i.M; };
    
    const T& lower() const { return m; };
    const T& upper() const { return M; };
    T& lower() { return m; };
    T& upper() { return M; };
    const T& inf() const { return lower(); };
    const T& sup() const { return upper(); };
    T& inf() { return m; };
    T& sup() { return M; };
    
    Interval& operator=( const T& x ) { m = M = x; return *this; };
    template<class X, int R>
    Interval& operator=( const Interval<X,R> & x );
    //    template<class X>
    //    Interval& operator+=( const X & x ) { *this = *this + x; return *this; };
    //    template<class X>
    //    Interval& operator-=( const X & x ) { *this = *this - x; return *this; };
    //    template<class X>
    //    Interval& operator*=( const X & x ) { *this = *this * x; return *this; };
    //    template<class X>
    //    Interval& operator/=( const X & x ) { *this = *this / x; return *this; };

    Interval& operator+=( const Interval & x ) { add(*this,x); return *this; };
    Interval& operator-=( const Interval & x ) { sub(*this,x); return *this; };
    Interval& operator*=( const Interval & x ) { mul(*this,x); return *this; };
    Interval& operator/=( const Interval & x ) { div(*this,x); return *this; };

  T width() const { return M-m; };
  T size() const { return width(); };
  operator T() const  {return (lower()+upper())/2; };
  T center() const { return (M+m)/T(2); };
  };


template<class OSTREAM, class T, int r> inline void
print(OSTREAM& os, const Interval<T,r>& a) {
  //os<<"[";print(os,a.lower());os<<",";mmx::print(os,a.upper());os<<"]";
}

template<class T, int r>
void hull( Interval<T,r>& v, const Interval<T,r>& a, const Interval<T,r>& b );
template<class T,  int r> inline
bool intersect( Interval<T,r>& result, 
        const Interval<T,r>& a, const Interval<T,r>& b );

} //end namespace mmx

#include <realroot/Interval_fcts.hpp>

namespace mmx {

namespace texp {

  template<class X, int r>  struct kernelof_< Interval<X,r> >  {
    typedef typename kernelof<X>::T T; };
}

namespace let {

  template<class X, class Y, int r0, int r1> inline void 
  assign( Interval<X,r0>& a, const Interval<Y,r1> & b )
  {
    assign(a.m,b.m);
    assign(a.M,b.M);
  };
};

template<class C,int R>
template<class X> inline
Interval<C,R> & Interval<C,R>::operator=( const texp::template_expression<X> & x ) {
  rnd_t rnd;        // rounding mode verification
  let::assign(*this,x);  // evaluation of template expression
  return *this;
}

template<class C,int R>
template<class X> inline
Interval<C,R>::Interval( const texp::template_expression<X> & x ) {
  rnd_t rnd;        // rounding mode verification
  std::cout<<"assign "<<x<<std::endl;
  let::assign(*this,x);  // evaluation of template expression
}

template<class C, int mode>
template<class X, int R>
Interval<C,mode>& Interval<C,mode>::operator=( const Interval<X,R> & x )
{
  let::assign(*this,x);
  return *this;
};

  namespace let {
    template<class A,class B> void assign(A& a, const B&b);
    
    template<class C, class T> inline void 
    assign( Interval<C>& a, const T & b ) {
      C mn,mx;
      {
    numerics::rounding<C> rnd( numerics::rounding<C>::rnd_up() );
         assign(mx,b);
      }
      
      {
    numerics::rounding<C> rhd( numerics::rounding<C>::rnd_dw() );
    assign(mn,b);
      }
      if( mx < mn ) std::swap(mn,mx);
      new (&a) Interval<C>(mn,mx);
    };
  }





} //end namespace mmx
    
#endif