interval.h

// ---------------------------------------------------------
//
//  interval.h
//  Tyson Brochu 2011
//
// ---------------------------------------------------------

#ifndef TUNICATE_INTERVAL_H
#define TUNICATE_INTERVAL_H

#include <cassert>
#include <fenv.h>
#include "interval_base.h"

class Interval;
typedef Interval IntervalType;

#ifndef DEBUG
#define VERIFY() (void)0;
#else
#define VERIFY() assert( -v[0] <= v[1] ); \
assert( v[0] == v[0] ); \
assert( v[1] == v[1] );
#endif

#ifndef DEBUG
#define CHECK_ROUNDING_MODE() (void)0;
#else
#define CHECK_ROUNDING_MODE() assert( fegetround( ) == FE_UPWARD );
#endif


// ----------------------------------------
//
// class Interval:
//
// Stores the interval [a,b] as [-a,b] internally.  With proper arithmetic operations, this
// allows us to use only FE_UPWARD and avoid switching rounding modes over and over.
//
// ----------------------------------------


class Interval : public IntervalBase
{

    static int s_previous_rounding_mode;

public:

   // Internal representation
   double v[2];

   Interval( double val );
   Interval( double left, double right );
   Interval();

   virtual ~Interval() {}

   virtual double stored_left() const;
   virtual double stored_right() const;

   virtual Interval& operator+=(const Interval &rhs);
   virtual Interval& operator-=(const Interval &rhs);
   virtual Interval& operator*=(const Interval &rhs);

   virtual Interval operator+(const Interval &other) const;
   virtual Interval operator-(const Interval &other) const;
   virtual Interval operator*(const Interval &other) const;
   // WARNING: Division bad and should be avoided in exact computations
   virtual Interval operator/(const Interval &other) const;

   virtual Interval operator-( ) const;

   static void begin_special_arithmetic();
   static void end_special_arithmetic();

};

inline void create_from_double( double a, Interval& out );


// ----------------------------------------

inline Interval::Interval( double val )
{
   v[0] = -val;
   v[1] = val;
   VERIFY();
}

// ----------------------------------------

inline Interval::Interval( double left, double right )
{
   assert( left <= right );
   v[0] = -left;
   v[1] = right;
   VERIFY();
}

// ----------------------------------------

inline Interval::Interval()
{
   v[0] = 0;
   v[1] = 0;
   VERIFY();
}

// ----------------------------------------

inline double Interval::stored_left() const
{
   return v[0];
}

// ----------------------------------------

inline double Interval::stored_right() const
{
   return v[1];
}


// ----------------------------------------

inline Interval& Interval::operator+=(const Interval &rhs)
{
   CHECK_ROUNDING_MODE();
   VERIFY();
   v[0] += rhs.v[0];
   v[1] += rhs.v[1];
   VERIFY();

   return *this;
}

// ----------------------------------------

inline Interval& Interval::operator-=( const Interval& rhs )
{
   CHECK_ROUNDING_MODE();
   v[0] += rhs.v[1];
   v[1] += rhs.v[0];
   VERIFY();
   return *this;
}

// ----------------------------------------

inline Interval& Interval::operator*=( const Interval& rhs )
{
   CHECK_ROUNDING_MODE();
   Interval p = (*this) * rhs;
   *this = p;
   return *this;
}

// ----------------------------------------

inline Interval Interval::operator+(const Interval &other) const
{
   CHECK_ROUNDING_MODE();
   double v0 = v[0] + other.v[0];
   double v1 = v[1] + other.v[1  ];
   return Interval(-v0, v1);
}

// ----------------------------------------

inline Interval Interval::operator-(const Interval &other) const
{
   CHECK_ROUNDING_MODE();
   double v0 = v[0] + other.v[1];
   double v1 = v[1] + other.v[0];
   return Interval(-v0, v1);
}

// ----------------------------------------

inline Interval Interval::operator*(const Interval &other) const
{
   CHECK_ROUNDING_MODE();

   double neg_a = v[0];
   double b = v[1];
   double neg_c = other.v[0];
   double d = other.v[1];

   Interval product;

   if ( b <= 0 )
   {
      if ( d <= 0 )
      {
         product.v[0] = -b * d;
         product.v[1] = neg_a * neg_c;
      }
      else if ( -neg_c <= 0 && 0 <= d )
      {
         product.v[0] = neg_a * d;
         product.v[1] = neg_a * neg_c;
      }
      else
      {
         product.v[0] = neg_a * d;
         product.v[1] = b * -neg_c;
      }
   }
   else if ( -neg_a <= 0 && 0 <= b )
   {
      if ( d <= 0 )
      {
         product.v[0] = b * neg_c;
         product.v[1] = neg_a * neg_c;
      }
      else if ( -neg_c <= 0 && 0 <= d )
      {
         product.v[0] = std::max( neg_a * d, b * neg_c );
         product.v[1] = std::max( neg_a * neg_c, b * d );
      }
      else
      {
         product.v[0] = neg_a * d;
         product.v[1] = b * d;
      }

   }
   else
   {
      if ( d <= 0 )
      {
         product.v[0] = b * neg_c;
         product.v[1] = -neg_a * d;
      }
      else if ( -neg_c <= 0 && 0 <= d )
      {
         product.v[0] = b * neg_c;
         product.v[1] = b * d;
      }
      else
      {
         product.v[0] = -neg_a * neg_c;
         product.v[1] = b * d;
      }
   }

   return product;

}

// WARNING: Division bad and should be avoided in exact computations
inline Interval Interval::operator/(const Interval &other) const
{
    CHECK_ROUNDING_MODE();
    assert(other.v[0] > 0 || other.v[1] < 0);
    return (*this) * Interval(1 / other.v[0], 1 / other.v[1]);
}

// ----------------------------------------

inline Interval Interval::operator-( ) const
{
   CHECK_ROUNDING_MODE();
   return Interval( -v[1], v[0] );
}

// ----------------------------------------

inline void Interval::begin_special_arithmetic()
{
   s_previous_rounding_mode = fegetround();
   fesetround( FE_UPWARD );
}

// ----------------------------------------

inline void Interval::end_special_arithmetic()
{
   fesetround( s_previous_rounding_mode );
}


// ----------------------------------------

inline void create_from_double( double a, Interval& out )
{
   out.v[0] = -a;
   out.v[1] = a;
}


#endif