interval_boost.h

Go to the documentation of this file.

// Interval Wrapper Class for BOOST -*- C++ -*-

// $Id: interval_boost.h 347 2008-02-06 15:20:01Z schilly $
// Copyright (C) 2001-2004 Hermann Schichl and Ferenc Domes
//
// This file is part of the COCONUT API.  This library
// is free software; you can redistribute it and/or modify it under the
// terms of the Library GNU General Public License as published by the
// Free Software Foundation; either version 2, or (at your option)
// any later version.

// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// Library GNU General Public License for more details.

// As a special exception, you may use this file as part of a free software
// library without restriction.  Specifically, if other files instantiate
// templates or use macros or inline functions from this file, or you compile
// this file and link it with other files to produce an executable, this
// file does not by itself cause the resulting executable to be covered by
// the Library GNU General Public License.  This exception does not however
// invalidate any other reasons why the executable file might be covered by
// the Library GNU General Public License.

/*
It must be compiled with -O0
*/

#ifndef _INTERVAL_BOOST_H_
#define _INTERVAL_BOOST_H_

#include <iostream>
#include <values.h>
#include <boost/numeric/interval.hpp>
#include <boost/numeric/interval/checking.hpp>
#include <boost/numeric/interval/compare.hpp>
#include <boost/numeric/interval/policies.hpp>

namespace coco
{

#define coconut_init_interval()

// checking template
template<class T>
struct checking_my
{ 
  static T pos_inf(){return COCO_INF;}
  static T neg_inf(){return -COCO_INF;}
  static T inf()
  {
    //assert(std::numeric_limits<T>::has_infinity);
    //return std::numeric_limits<T>::infinity();
    return COCO_INF;
  }
  static T nan()
  {
    assert(std::numeric_limits<T>::has_quiet_NaN);
    return std::numeric_limits<T>::quiet_NaN();
  }
  static bool is_nan(const T& x)
  {
    return std::numeric_limits<T>::has_quiet_NaN && (x != x);
  }
  static T empty_inf()
  {
    return (std::numeric_limits<T>::has_quiet_NaN ?
            std::numeric_limits<T>::quiet_NaN() : static_cast<T>(1));
  }
  static T empty_upper()
  {
    return (std::numeric_limits<T>::has_quiet_NaN ?
            std::numeric_limits<T>::quiet_NaN() : static_cast<T>(0));
  }
  static T empty_lower()
  {
    return (std::numeric_limits<T>::has_quiet_NaN ?
            std::numeric_limits<T>::quiet_NaN() : static_cast<T>(0));
  }
  static bool is_empty(const T& l, const T& u)
  {
    return !(l <= u); // safety for partial orders
  }
};

template<class T>
struct my_rounded_math:
  boost::numeric::interval_lib::save_state_nothing<
    boost::numeric::interval_lib::rounded_transc_opp<T>
    >
{};

//policies type
typedef boost::numeric::interval_lib::policies<
        my_rounded_math<double>,
        checking_my<double>
>my_policies;

//data type
class interval;

struct interval_st
{
  double l, u;
  interval_st& operator=(const interval& __i);
};

//CLASS BEGIN
class interval
{
private:
//interval type
  typedef boost::numeric::interval<double, my_policies> __Base;

private:
  __Base I;

//unprotect the interval type
// save and initialize the rounding mode
//typename bointerval_prot::traits_type::rounding rnd;
// define the unprotected version of the interval type
//typedef boost::numeric::interval_lib::unprotect<bointerval_prot>::type bointerval;

  static int precision_val;
  
private:
  static interval __construct(const __Base& x)
        { interval __x; __x.I = x; return __x; }

public:
  interval(double lo, double up) : I(lo, up) {}
  interval(double d = 0) : I(d) {}
  interval(int d) : I((double) d) {}
  interval(unsigned d) : I((double) d) {}
  interval(long d) : I((double) d) {}
  interval(unsigned long d) : I((double) d) {}
  interval(const interval& x) : I(x.base()) {}
  interval(const interval_st& x) : I(x.l, x.u) {}

  //i.set(a,b) sets i to [a, b]
  void set(double lo, double up)
  {
    *this = interval(lo, up);
  }
  
  //EXTRA BEGIN
  double inf() const {return I.lower();}
  double sup() const {return I.upper();}
  __Base base() const {return I;}
  static interval EMPTY() {return __construct(__Base::empty());}
  //END EXTRA
  
  bool empty() const {return boost::numeric::empty(I);}
  bool is_empty() const {return  boost::numeric::empty(I);}
  bool is_thin() const {return singleton(I);}
  bool is_unbounded_below() const {return inf() == -COCO_INF;}
  bool is_unbounded_above() const {return sup() == COCO_INF;}
  bool is_entire() const {return inf() == -COCO_INF && sup() == COCO_INF;}
  bool is_bounded() const
        {return !is_unbounded_below() && !is_unbounded_above();}
  //EXTRA BEGIN
  bool subset (const interval& x) const 
        {return boost::numeric::subset(I, x.base());}
  bool superset (const interval& x) const 
        {return boost::numeric::subset(x.base(),I);}
  bool proper_subset (const interval& x) const 
        {return boost::numeric::proper_subset(I, x.base());}
  bool proper_superset (const interval& x) const 
        {return boost::numeric::proper_subset(x.base(),I);}
  bool disjoint (const interval& x) const 
        {return !boost::numeric::overlap(I, x.base());}
  //END EXTRA
  double rel_width() const;

  // y.intersect(x) intersects y with x
  interval& intersect(const interval& x)
  {
    I=boost::numeric::intersect(I,x.base());
    return *this;
  }

  // y.hull(x) builds the hull of y with x
  interval& hull(const interval& x)
  {
    I=boost::numeric::hull(I,x.base());
    return *this;
  }

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

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

  void setpair(double& lo, double& up) const {lo = inf(); up = sup();}
  void get_bounds(double& a, double& b) const {a = inf(); b = sup();}
  void set_bounds(double a, double b) {I.assign(a, b);}
  void set_lb(double d) {I.assign(d, I.upper());}
  void set_ub(double d) {I.assign(I.lower(), d);}
  
  interval& operator=(const interval& x){I=__Base(x.base()); return *this;}
  interval& operator=(double d) {I.assign(d,d); return *this;}
  interval& operator=(int d) {I.assign(d,d); return *this;}
  interval& operator=(unsigned d) {I.assign(d,d); return *this;}
  interval& operator=(long d) {I.assign(d,d); return *this;}
  interval& operator=(unsigned long d) {I.assign(d,d); return *this;}
  //interval& operator=(const __Base& _x)
  //        {I.assign(_x.lower(),_x.upper()); return *this;}

  interval operator-() const {return __construct(-I);}

  interval& ipow(int n);
  interval& imin(const interval& x);
  interval& imax(const interval& x);

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

  //--------------------------------------------------------------------------
  // FUNCTIONS THAT NEED TO BE PUT IN CPP FILE
  //--------------------------------------------------------------------------

  
  //unary +=
  
  interval& operator+=(const interval& a);

  
  //unary += for number
  
  interval& operator+=(double a);

  
  //unary -=
  
  interval& operator-=(const interval& a);

  //unary -= for number
  
  interval& operator-=(double a);

  //unary operator*=

  interval& operator*=(const interval& a);

  //unary operator*= for point
  
  interval& operator*=(double a);

  //unary operator/=
  
  interval& operator/=(const interval& a);

  //unary operator/= for point

  interval& operator/=(double a);

  double mid() const;

  /*Returns an upper bound for the diameter (width) of this interval, i.e.

  a.diam() == a.sup()-a.inf()

  Special cases in the extended system:
  <UL>
  <LI> a.diam() == NaN for a == [ EMPTY ]
  <LI> a.diam() == +INF for any infinite interval
  </UL>
  */
  
  double diam() const;
  double width() const;

  double relDiam() const;

  double rad() const;

  double mig() const;

  double mag() const;

  double dist(const interval& x) const;

  double safeguarded_mid() const;

  static int const & precision();
  static int precision(int const & p);

  friend std::ostream &operator << (std::ostream &s, const interval &a);
  friend interval operator+(const interval& a, const interval& b);
  friend interval operator+(const interval& a, double b);
  friend interval operator+(double b, const interval& a);
  friend interval operator-(const interval& a, const interval& b);
  friend interval operator-(const interval& a, double b);
  friend interval operator-(double b, const interval& a);
  friend interval cancel(const interval& a, const interval& b);
  friend interval operator*(const interval& a, const interval& b);
  friend interval operator*(const interval& a, double b);
  friend interval operator*(double b, const interval& a);
  friend interval operator/(const interval& a, const interval& b);
  friend interval operator/(const interval& a, double b);
  friend interval operator/(double b, const interval& a);
  friend interval acos(const interval& x);
  friend interval abs(const interval& x);
  friend interval acosh(const interval& x);
  friend interval acot(const interval& x);
  friend interval acoth(const interval& x);
  friend interval asin(const interval& x);
  friend interval asinh(const interval& x);
  friend interval atan(const interval& x);
  friend interval atanh(const interval& x);
  friend interval cos(const interval& x);
  friend interval cosh(const interval& x);
  friend interval cot(const interval& x);
  friend interval coth(const interval& x);
  friend interval exp(const interval& x);
  friend interval exp10(const interval& x);
  friend interval exp2(const interval& x);
  friend interval expm1(const interval& x);
  friend interval log(const interval& x);
  friend interval log10(const interval& x);
  friend interval log1p(const interval& x);
  friend interval log2(const interval& x);
  friend interval power(const interval& x, int n);
  friend interval pow(const interval& x, const interval& y);
  friend interval sin(const interval& x);
  friend interval sinh(const interval& x);
  friend interval sqr(const interval& x);
  friend interval sqrt(const interval& x);
  friend interval tan(const interval& x);
  friend interval tanh(const interval& x);
  friend interval imax(const interval& x, const interval& y);
  friend interval imin(const interval& x, const interval& y);
  friend interval division_part1(const interval& x, const interval& y,
        bool& b);
  friend interval division_part2(const interval& x, const interval& y,
        bool b);
};
// end of interval class


//------------------------------------------------------------------------------
// INLINE FUNCTIONS
//------------------------------------------------------------------------------

inline double safeguarded_mid(const interval& __i)
        {return __i.safeguarded_mid();}
interval ipow(const interval& x, int n);

interval gauss(const interval& x);
interval atan2(const interval& y, const interval& x);
double absmin(const interval& __i);
double gainfactor(const interval& _old, const interval& _new);

template <class _TC>
bool operator==(const interval& __i, const _TC& __d);

template <class _TC>
bool operator!=(const interval& __i, const _TC& __d);



//------------------------------------------------------------------------------
// FUNCTIONS TO BE HIDDEN FROM USERS TO AVOID THE OPTIMIZATION OF COMPILER
//------------------------------------------------------------------------------

//EXTRA COMPARSION OPERATORS
inline bool operator==(const interval& a, const interval& b)
        {return a.sup()==b.sup() && a.inf()==b.inf();}
inline bool operator==(const interval& a, double b)
        {return a.sup()==b && a.inf()==b;}
inline bool operator!=(const interval& a, const interval& b)
        {return a.sup()!=b.sup() || a.inf()!=b.inf();}
inline bool operator!=(const interval& a, double b)
        {return a.sup()!=b || a.inf() != b;}
inline bool operator<(const interval& a, const interval& b){return a.sup() < b.inf();}
inline bool operator<(const interval& a, double b){return a.sup() < b;}
inline bool operator<(double a, const interval& b){return a < b.inf();}
//END EXTRA
interval operator+(const interval& a, const interval& b);
interval operator+(const interval& a, double b);
interval operator+(double b, const interval& a);
interval operator-(const interval& a, const interval& b);
interval operator-(const interval& a, double b);
interval operator-(double b, const interval& a);
interval cancel(const interval& a, const interval& b);
interval operator*(const interval& a, const interval& b);
interval operator*(const interval& a, double b);
interval operator*(double b, const interval& a);
interval operator/(const interval& a, const interval& b);
interval operator/(const interval& a, double b);
interval operator/(double b, const interval& a);

//OUTPUT
std::ostream &operator << (std::ostream &s, const interval &a);

double mid(const interval&);
double diam(const interval&);
double width(const interval&);
double relDiam(const interval&);
double rad(const interval&);
double mig(const interval&);
double mag(const interval&);
double dist(const interval& x, const interval& y);

interval round_to_integer(const interval& x);

interval acos(const interval& x);

interval abs(const interval& x);

interval acosh(const interval& x);

interval acot(const interval& x);

interval acoth(const interval& x);

interval asin(const interval& x);

interval asinh(const interval& x);

interval atan(const interval& x);

interval atanh(const interval& x);

interval cos(const interval& x);

interval cosh(const interval& x);

interval cot(const interval& x);

interval coth(const interval& x);

interval exp(const interval& x);

interval exp10(const interval& x);

interval exp2(const interval& x);

interval expm1(const interval& x);

interval log(const interval& x);

interval log10(const interval& x);

interval log1p(const interval& x);

interval log2(const interval& x);

interval power(const interval& x, int n);

interval pow(const interval& x, const interval& y);

interval sin(const interval& x);

interval sinh(const interval& x);

interval sqr(const interval& x);

interval sqrt(const interval& x);

interval tan(const interval& x);

interval tanh(const interval& x);

interval imax(const interval& x, const interval& y);

interval imin(const interval& x, const interval& y);


//------------------------------------------------------------------------------
// INLINE IMPLEMENTAION
//------------------------------------------------------------------------------

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


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


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


inline
interval atan2(const interval& y, const interval& x)
{
  return atan(y/x);
}


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

interval division_part1(const interval& x, const interval& y, bool& b);

interval division_part2(const interval& x, const interval& y, bool b = true);

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() != COCO_INF && _old.sup() == COCO_INF) ||
     (_new.inf() != -COCO_INF && _old.inf() == -COCO_INF))
    return 0.;
  return 1.;
}


template <class _TC>
inline
bool operator==(const interval& __i, const _TC& __d)
{
  return boost::numeric::interval_lib::compare::set::operator==(__i.base(), interval(__d).base());
}


template <class _TC>
inline
bool operator!=(const interval& __i, const _TC& __d)
{
  return boost::numeric::interval_lib::compare::set::operator!=(__i.base(), interval(__d).base());
}


// special additional functions
inline
double interval::project(double __d) const
{
  if(__d < inf())
    return inf();
  else if(__d > sup())
    return sup();
  else
    return __d;
}


} // end of namespace coco


#endif // _INTERVAL_BOOST_H_

---