---

interval.h

Go to the documentation of this file.

#ifndef _INTERVAL_H_
#define _INTERVAL_H_

#define FILIB_USE_MACROS 0
#define FILIB_EXTENDED 1

#include <iostream>
#if FILIB_USE_MACROS
#include <Interval.h>
#else
#include <interval/interval.hpp>
#endif

#if defined(INFINITY)
#undef INFINITY
#endif
#define INFINITY filib::fp_traits_base<double>::infinity()

#if defined(L_PI)
#undef L_PI
#endif
#define L_PI filib::fp_traits_base<double>::l_pi()

#if defined(U_PI)
#undef U_PI
#endif
#define U_PI filib::fp_traits_base<double>::u_pi()

//#define min(a,b) ((a)<(b)?(a):(b))
//#define max(a,b) ((a)>(b)?(a):(b))
#define sgn(x)   ((x)>0?1.:(x)<0?-1.:0.)
#define is_integer(x) (rint((x))==(x)&&rint((x))<MAXINT&&rint((x))>-MAXINT)
#define get_integer(x) ((int)rint((x)))

//data type
class interval;

struct interval_st
{
  double l,u;

  interval_st& operator=(const interval& __i);
};

#if FILIB_USE_MACROS
class interval : public Interval
#else
class interval : public filib::interval<>
#endif
{
private:
#if FILIB_USE_MACROS
  typedef Interval __Base;
#else
  typedef filib::interval<> __Base;
#endif

public:
  interval(const double& lo, const double& up) : __Base(lo,up) {}

  interval(const double& __d = 0) : __Base(__d) {}
  
  interval(const int& __d) : __Base(__d+0.) {}
  interval(const unsigned int& __d) : __Base(__d+0.) {}
  interval(const long int& __d) : __Base(__d+0.) {}
  interval(const unsigned long int& __d) : __Base(__d+0.) {}

  interval(const __Base& __i) : __Base(__i) {}

  interval(const interval& __i) : __Base(__i) {}

  interval(const interval_st& __i) : __Base(__i.l, __i.u) {}
  
  //i.set(a,b) sets i to [a,b]
  void set(double lo, double up)
  {
    *this = interval(lo,up);
  }

  bool empty() const { return __Base::isEmpty(); }
  bool is_thin() const { return __Base::isPoint(); }
  bool is_unbounded_below() const { return inf() == -INFINITY; }
  bool is_unbounded_above() const { return sup() == INFINITY; }
  bool is_entire() const { return inf() == -INFINITY && sup() == INFINITY; }
  bool is_bounded() const
        { return !is_unbounded_below() && !is_unbounded_above(); }

  double rel_width() const { return __Base::relDiam(); }

  // i.intersect(j) intersects i with j
  interval& intersect(const interval& __i)
  {
    interval help(__Base::intersect(__i));
    return *this = help;
  }

  interval& hulldiff(const interval& __i)
  {
    if(subset(__i))
      *this = EMPTY();
    else if(!superset(__i))
    {
      if(__i.contains(inf()))
        *this = interval(__i.sup(), sup());
      else
        *this = interval(inf(), __i.inf());
    }
    return *this;
  }

  bool contains(const interval& __i) const
  {
    if(inf() <= __i.inf() && sup() >= __i.sup())
      return true;
    else
      return false;
  }

  void setpair(double& __l, double& __u) const { __l = inf(); __u = sup(); }

  interval& operator=(const double& __d) { return *this = interval(__d); }
  interval& operator=(const int& __d) { return *this = interval(__d); }
  interval& operator=(const unsigned int& __d) { return *this = interval(__d); }
  interval& operator=(const long int& __d) { return *this = interval(__d); }
  interval& operator=(const unsigned long int& __d)
        { return *this = interval(__d); }

  interval& ipow(int n) { *this = power(*this, n); return *this; }

  interval& intersect_div(const interval& __i, const interval& __j);
  interval& intersect_power(const interval& __i, const int& __n);
  interval& intersect_invsqr_wc(const interval& __i, const double& __l,
                                const double& __d);
  interval& intersect_invpower_wc(const interval& __i, const double& __l,
                               const int& __n);
  interval& intersect_invsin_wc(const interval& __i, const double& __l,
                                const double& __d);
  interval& intersect_invcos_wc(const interval& __i, const double& __l,
                                const double& __d);
  interval& intersect_invgauss_wc(const interval& __i, const double& __l,
                                  const double& __m, const double& __s);

  double project(const double& __d);
};
// end of interval class

inline interval_st& interval_st::operator=(const interval& __i)
{ 
  l=__i.inf();
  u=__i.sup();
  return *this;
}

// special 'splitting' methods
inline interval& interval::intersect_div(const interval& __i,
                                         const interval& __j)
{
  interval help;
  
  if(__j.inf() < 0 && __j.sup() > 0 && !__i.contains(0))
  {
    interval up, lo;
    
    up = __i/interval(0,__j.sup());
    lo = __i/interval(__j.inf(),0);
    up.intersect(*this);
    lo.intersect(*this);
    if(up.empty()) *this = lo;
    else if(lo.empty()) *this = up;
  }
  else
  {
    help = __i/__j;
    this->intersect(help);
  }
  return *this;
}

inline interval& interval::intersect_power(const interval& __i,
                                           const int& __n)
{
  if(__n < 0 && __i.inf() < 0 && __i.sup() > 0)
  {
    if((-__n)&1) // __n odd
    {
      interval up, lo;
      
      up = power(interval(0,__i.sup()),__n);
      lo = power(interval(__i.inf(),0),__n);
      up.intersect(*this);
      lo.intersect(*this);
      if(up.empty()) *this = lo;
      else if(lo.empty()) *this = up;
    }
    else // __n even
    {
      interval res;

      res = 1./power(__i, -__n);
      this->intersect(res);
    }
  }
  else
    this->intersect(power(__i,__n));
  return *this;
}

inline interval& interval::intersect_invsqr_wc(const interval& __i,
                        const double& __l, const double& __d)
{
  interval up, lo;

  up = (sqrt(__i) - __d)/__l;
  lo = (-sqrt(__i) - __d)/__l;
  up.intersect(*this);
  lo.intersect(*this);
  if(up.empty())
    *this = lo;
  else if(lo.empty())
    *this = up;
  else
    this->intersect(interval(lo.inf(), up.sup()));
  return *this;
}

inline interval& interval::intersect_invpower_wc(const interval& __i,
                        const double& __l, const int& __n)
{
  if(__n == 0)
  {
    if(!__i.contains(1.))
      *this = EMPTY();
  }
  else if(__n > 0)
  {
    if(__n & 1) // n is odd, hence power bijective R^*<->R^*
    {
      interval __exp = interval(1.)/(1.*__n);
      interval __res;

      if(__i.inf() < 0 && __i.sup() > 0)
      {
        __res = -pow(interval(0,-__i.inf()),__exp);
        __res.hull(pow(interval(0,__i.sup()),__exp));
      }
      else if(__i.inf() < 0)
        __res = -pow(-__i,__exp);
      else
        __res = pow(__i,__exp);
      this->intersect(__res/__l);
    }
    else
    {
      interval __exp = interval(1.)/(1.*__n);
      interval __hi, __lo;

      if(__i.sup() < 0)
      {
        *this = EMPTY();
        return *this;
      }
      else if(__i.inf() < 0)
        __hi = pow(interval(0,__i.sup()),__exp)/__l;
      else
        __hi = pow(__i,__exp)/__l;
      __lo = -__hi;
      __hi.intersect(*this);
      __lo.intersect(*this);
      if(__lo.empty())
        *this = __hi;
      else if(__hi.empty())
        *this = __lo;
      else
        *this = interval(__lo.inf(), __hi.sup());
    }
  }
  else // __n < 0
  {
    if((-__n)&1) // __n odd
    {
      interval __exp = interval(1.)/(-1.*__n);
      interval __res;

      if(__i.inf() < 0 && __i.sup() > 0)
      {
        interval __up;
        
        __res = -pow(interval(0,-__i.inf()),__exp);
        __res.hull(pow(interval(0,__i.sup()),__exp));
        __res *= __l;
        __up = 1./interval(0,__res.sup());
        __res = 1./interval(__res.inf(),0);
        __up.intersect(*this);
        __res.intersect(*this);
        if(__up.empty()) *this = __res;
        else if(__res.empty()) *this = __up;
        // else *this is unchanged, because hull(__res,__up) == ENTIRE()
      }
      else if(__i.inf() < 0)
      { // the argument interval is cle 0
        __res = -pow(-__i,__exp)*__l;
        this->intersect(1./__res);
      }
      else
      { // the argument interval is cge 0
        __res = pow(__i,__exp)*__l;
        this->intersect(1./__res);
      }
    }
    else // __n is even and < 0
    {
      interval __exp = interval(1.)/(1.*__n);
      interval __up, __lo;

      if(__i.sup() < 0)
        *this = EMPTY();
      else if(__i.inf() < 0)
        __lo = interval(0,__i.sup());
      else
        __lo = __i;

      __up = pow(__lo,__exp)/__l;
      __lo = -__up;
      __lo.intersect(*this);
      __up.intersect(*this);
      if(__lo.empty())
        *this = __up;
      else if(__up.empty())
        *this = __lo;
      else
        *this = interval(__lo.inf(), __up.sup());
    }
  }
  return *this;
}

inline interval& interval::intersect_invsin_wc(const interval& __i,
                        const double& __l, const double& __d)
{
  std::cerr << "interval intersect_invsin_wc is not yet implemented!" <<
    std::endl;
  exit(-1);
  return *this;
}

inline interval& interval::intersect_invcos_wc(const interval& __i,
                        const double& __l, const double& __d)
{
  std::cerr << "interval intersect_invcos_wc is not yet implemented!" <<
    std::endl;
  exit(-1);
  return *this;
}

inline interval& interval::intersect_invgauss_wc(const interval& __i,
                const double& __l, const double& __m, const double& __s)
{
  std::cerr << "interval intersect_invgauss_wc is not yet implemented!" <<
    std::endl;
  exit(-1);
  return *this;
}

// elementary functions for intervals
inline interval ipow(const interval& __i, int __n)
{
  return power(__i,__n);
}

inline interval gauss(const interval& __i)
{
  return exp(-sqr(__i));
}

inline interval atan2(const interval& __i, const interval& __j)
{
  std::cerr << "interval atan2 is not yet implemented!" << std::endl;
  exit(-1);
}

inline double safeguarded_mid(const interval& __i)
{
  if(__i.inf() < 0 && __i.sup() > 0)
    return 0.;
  else if(__i.inf() == -INFINITY)
    return 2.*__i.sup();
  else if(__i.sup() == INFINITY)
    return 2.*__i.inf();
  else if(__i.inf() == 0 || __i.sup() == 0)
    return mid(__i);
  else
  {
    double r = sgn(__i.inf())*sqrt(fabs(__i.sup()))*sqrt(fabs(__i.inf()));
    if(__i.contains(r))
      return r;
    else
      return mid(__i);
  }
}

inline double absmin(const interval& __i)
{
  if(__i.contains(0.))
    return 0.;
  else if(__i.inf() > 0)
    return __i.inf();
  else
    return __i.sup();
}

template <class _C>
inline bool operator==(const interval& __i, const _C& __d)
{
  return filib::operator==(__i, interval(__d));
}

template <class _C>
inline bool operator!=(const interval& __i, const _C& __d)
{
  return filib::operator!=(__i, interval(__d));
}

inline double gainfactor(const interval& _old, const interval& _new)
{
  if(_old.is_thin())
    return 1.;
  if(_old.is_bounded() && _new.is_bounded())
    return _new.width()/_old.width();
  if((_new.sup() != INFINITY && _old.sup() == INFINITY) ||
     (_new.inf() != -INFINITY && _old.inf() == -INFINITY))
    return 0.;
  return 1.;
}

#if 0
inline interval operator+(int n, const interval& i) { return interval(n)+i; }
inline interval operator+(const interval& i, int n) { return interval(n)+i; }
inline interval operator-(int n, const interval& i) { return interval(n)-i; }
inline interval operator-(const interval& i, int n) { return interval(n)-i; }
inline interval operator*(int n, const interval& i) { return interval(n)*i; }
inline interval operator*(const interval& i, int n) { return interval(n)*i; }
inline interval operator/(int n, const interval& i) { return interval(n)/i; }
inline interval operator/(const interval& i, int n) { return interval(n)/i; }
#endif

// special additional functions
#if 0
double project(const double& __d) const
{
  if(__d < inf())
    return inf();
  else if(__d > sup())
    return sup();
  else
    return __d;
}
#endif
#endif // _INTERVAL_H_

---