/Users/mourrain/Devel/mmx/numerix/include/numerix/interval.hpp

Go to the documentation of this file.

/******************************************************************************
* MODULE     : interval.hpp
* DESCRIPTION: Intervals
* COPYRIGHT  : (C) 2006  Joris van der Hoeven
*******************************************************************************
* This software falls under the GNU general public license and comes WITHOUT
* ANY WARRANTY WHATSOEVER. See the file $TEXMACS_PATH/LICENSE for more details.
* If you don't have this file, write to the Free Software Foundation, Inc.,
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
******************************************************************************/

#ifndef __MMX_INTERVAL_HPP
#define __MMX_INTERVAL_HPP
#include <basix/function.hpp>
#include <numerix/rounded.hpp>
namespace mmx {
#define TMPL_DEF template<typename C,typename V=std_interval>
#define TMPL template<typename C,typename V>
#define Interval interval<C,V>
#define Up V::template rnd_helper<C>::UV
#define Down V::template rnd_helper<C>::DV
TMPL class interval;
TMPL inline C lower (const Interval& z);
TMPL inline C upper (const Interval& z);

/******************************************************************************
* Interval variants
******************************************************************************/

struct std_interval {
  template<typename C>
  struct rnd_helper {
    typedef Round_up(C) UV;
    typedef Round_down(C) DV;
  };
};

/******************************************************************************
* The interval type
******************************************************************************/

TMPL_DEF
class interval {
MMX_ALLOCATORS
private:
  C l;
  C r;

public:      
  typedef C boundary_type;

  inline interval (): l (0), r (0) {}
  template<typename T> inline interval (const T& x):
    l (Down::template as<C> (x)),
    r (Up::template as<C> (x)) {}
  template<typename T> inline interval (const interval<T,V>& x):
    l (Down::template as<C> (lower (x))),
    r (Up::template as<C> (upper (x))) {}
  template<typename LT,typename UT> inline
  interval (const LT& l2, const UT& u2):
    l (Down::template as<C> (l2)),
    r (Up::template as<C> (u2)) {}
  inline interval (const C& x): l (x), r (x) {}
  inline interval (const C& l2, const C& r2): l (l2), r (r2) {}
  friend C lower LESSGTR (const Interval& x);
  friend C upper LESSGTR (const Interval& x);

  static inline Interval nan ();
  static inline Interval maximal ();
  static inline Interval minimal ();
  static inline Interval infinity ();
  static inline Interval fuzzy ();
  static inline Interval log2 ();
  static inline Interval pi ();
  static inline Interval euler ();

  inline Interval& operator += (const Interval& x);
  inline Interval& operator -= (const Interval& x);
  inline Interval& operator <<= (const xint& shift);
  inline Interval& operator >>= (const xint& shift);
};
DEFINE_UNARY_FORMAT_1 (interval)

UNARY_RETURN_TYPE(TMPL,center_op,Interval,C);
UNARY_RETURN_TYPE(TMPL,radius_op,Interval,C);
TMPL struct rounding_helper<Interval >: public rounding_helper<C> {};

/******************************************************************************
* Standard constructors and accessors
******************************************************************************/

template<typename C, typename V, typename C2, typename R2>
struct make_ball_helper<Interval,C2,R2> {
  static inline Interval val (const C2& c, const R2& r) {
    typedef Round_up(C2) Rnd_Up;
    typedef Round_down(C2) Rnd_Down;
    return Interval (Rnd_Down::template sub_as<C2> (c, r),
                     Rnd_Up::template add_as<C2> (c, r));
  }
};

template<typename C, typename V, typename C2>
struct make_interval_helper<Interval,C2> {
  static inline Interval val (const C2& l, const C2& r) {
    return Interval (l, r); }
};

TMPL inline C lower (const Interval& x) { return x.l; }
TMPL inline C upper (const Interval& x) { return x.r; }
TMPL inline C center (const Interval& x) {
  return decexp2 (lower (x) + upper (x)); }
TMPL inline C radius (const Interval& x) {
  if (is_infinite (x)) return 0;
  return Up::add (decexp2 (Up::sub (upper (x), lower (x))),
                  multiplicative_error (center (x))); }
TMPL inline format<C> CF (const Interval& x) {
  return get_format (lower (x)); }

TMPL inline Interval abs (const Interval& x);
TMPL inline C abs_down (const Interval& x) { return lower (abs (x)); }
TMPL inline C abs_up (const Interval& x) { return upper (abs (x)); }
TMPL inline C bnd_down (const Interval& x) { return lower (x); }
TMPL inline C bnd_up (const Interval& x) { return upper (x); }

/******************************************************************************
* Basic routines
******************************************************************************/

TMPL inline Interval copy (const Interval& x) {
  return Interval (copy (lower (x)), copy (upper (x))); }
TMPL inline Interval duplicate (const Interval& x) {
  return Interval (duplicate (lower (x)), duplicate (upper (x))); }

TMPL inline bool is_exact_zero (const Interval& x) {
  return is_exact_zero (lower (x)) && is_exact_zero (upper (x)); }

template<typename Op, typename C, typename V> nat
unary_hash (const Interval& x) {
  nat h= Op::op (lower (x));
  return (h<<3) ^ (h<<15) ^ (h>>29) ^ Op::op (upper (x));
}

template<typename Op, typename C, typename V> bool
binary_test (const Interval& x1, const Interval& x2) {
  return Op::op (lower (x1), lower (x2)) && Op::op (upper (x1), upper (x2));
}

template<typename S1, typename D> interval<D>
map (const function_1<D,Argument(S1) >& fun, const interval<S1>& x) {
  return interval<D> (fun (lower (x)), fun (upper (x)));
}

template<typename S1, typename D> interval<D>
map (D (*fun) (const S1&), const interval<S1>& x) {
  return interval<D> (fun (lower (x)), fun (upper (x)));
}

EXACT_IDENTITY_OP_SUGAR(TMPL,Interval)
HARD_IDENTITY_OP_SUGAR(TMPL,Interval)
TMPL nat hash (const Interval& x) { return 12721; }

TMPL syntactic
flatten (const Interval& x) {
  if (is_nan (x)) return Nan (syntactic);
  if (lower (x) == 0 && upper (x) == 0) return 0;
  return flatten_range (lower (x), upper (x));
  // return flatten ("interval", lower (x), upper (x));
}

TMPL syntactic
flatten (const Interval& x, xnat dd) {
  if (is_nan (x)) return Nan (syntactic);
  if (lower (x) == 0 && upper (x) == 0) return 0;
  return flatten_range (lower (x), upper (x), dd);
}

TMPL
struct binary_helper<Interval >: public void_binary_helper<Interval > {
  static inline string short_type_name () {
    return "I" * Short_type_name (C); }
  static inline generic full_type_name () {
    return gen ("Interval", Full_type_name (C)); }
  static inline generic disassemble (const Interval& x) {
    return gen_vec (as<generic> (lower (x)), as<generic> (upper (x))); }
  static inline Interval assemble (const generic& v) {
    return Interval (as<C> (vector_access (v, 0)),
                     as<C> (vector_access (v, 1))); }
  static inline void write (const port& out, const Interval& x) {
    binary_write<C> (out, lower (x));
    binary_write<C> (out, upper (x)); }
  static inline Interval read (const port& in) {
    C l= binary_read<C> (in);
    C r= binary_read<C> (in);
    return Interval (l, r); }
};

/******************************************************************************
* Comparisons
******************************************************************************/

TMPL inline bool
intersect (Interval& r, const Interval &a, const Interval &b) {
  if ((upper (a) < lower (b)) || (lower (a) > upper (b))) return false;
  r = Interval (max (lower (a), lower (b)), min (upper (a), upper (b)));
  return true;
}

TMPL inline bool
is_zero (const Interval& x) {
  return is_finite (x) && sign (lower (x)) * sign (upper (x)) <= 0;
}

TMPL inline bool
is_non_zero (const Interval& x) {
  return !is_finite (x) || sign (lower (x)) * sign (upper (x)) > 0;
}

TMPL inline bool
operator == (const Interval& x1, const Interval& x2) {
  if (is_nan (x1) || is_nan (x2)) return is_nan (x1) && is_nan (x2);
  return lower (x1) <= upper (x2) && lower (x2) <= upper (x1);
}

TMPL inline bool
operator != (const Interval& x1, const Interval& x2) {
  if (is_nan (x1) || is_nan (x2)) return !(is_nan (x1) && is_nan (x2));
  return lower (x1) > upper (x2) || lower (x2) > upper (x1);
}

TMPL inline bool
operator <= (const Interval& x1, const Interval& x2) {
  if (is_nan (x1) || is_nan (x2)) return false;
  return lower (x1) <= upper (x2);
}

TMPL inline bool
operator < (const Interval& x1, const Interval& x2) {
  if (is_nan (x1) || is_nan (x2)) return false;
  return lower (x2) > upper (x1);
}

TMPL inline bool
operator >= (const Interval& x1, const Interval& x2) {
  if (is_nan (x1) || is_nan (x2)) return false;
  return lower (x2) <= upper (x1);
}

TMPL inline bool
operator > (const Interval& x1, const Interval& x2) {
  if (is_nan (x1) || is_nan (x2)) return false;
  return lower (x1) > upper (x2);
}

EQUAL_INT_SUGAR(TMPL,Interval);
COMPARE_INT_SUGAR(TMPL,Interval);

/******************************************************************************
* Other order related functions
******************************************************************************/

TMPL inline Interval
abs (const Interval& x) {
  if (is_nan (x)) return x;
  C r= max (abs (lower (x)), abs (upper (x)));
  if (sign (lower (x)) * sign (upper (x)) <= 0) return Interval (C(0), r);
  return Interval (min (abs (lower (x)), abs (upper (x))), r);
}

TMPL Interval
min (const Interval& x1, const Interval& x2) {
  return Interval (min (lower (x1), lower (x2)), min (upper (x1), upper (x2)));
}

TMPL Interval
max (const Interval& x1, const Interval& x2) {
  return Interval (max (lower (x1), lower (x2)), max (upper (x1), upper (x2)));
}

TMPL Interval
floor (const Interval& x) {
  return Interval (Down::floor (lower (x)), Up::floor (upper (x)));
}

TMPL Interval
trunc (const Interval& x) {
  return Interval (Down::trunc (lower (x)), Up::trunc (upper (x)));
}

TMPL Interval
ceil (const Interval& x) {
  return Interval (Down::ceil (lower (x)), Up::ceil (upper (x)));
}

TMPL Interval
round (const Interval& x) {
  return Interval (Down::round (lower (x)), Up::round (upper (x)));
}

/******************************************************************************
* Constants
******************************************************************************/

TMPL inline void set_nan (Interval& x) {
  x= Interval (Nan (C), Nan (C)); }
TMPL inline void set_maximal (Interval& x) {
  x= Interval (Maximal (C), Maximal (C)); }
TMPL inline void set_minimal (Interval& x) {
  x= Interval (Minimal (C), Minimal (C)); }
TMPL inline void set_infinity (Interval& x) {
  x= Interval (Infinity (C), Infinity (C)); }
TMPL inline void set_fuzz (Interval& x) {
  x= Interval (Minimal (C), Maximal (C)); }
TMPL inline void set_smallest (Interval& x) {
  x= Interval (Smallest (C), Smallest (C)); }
TMPL inline void set_largest (Interval& x) {
  x= Interval (Largest (C), Largest (C)); }
TMPL inline void set_accuracy (Interval& x) {
  x= Interval (Accuracy (C), Accuracy (C)); }
TMPL inline void set_log2 (Interval& x) {
  x= Interval (Down::template log2<C>(), Up::template log2<C>()); }
TMPL inline void set_pi (Interval& x) {
  x= Interval (Down::template pi<C>(), Up::template pi<C>()); }
TMPL inline void set_euler (Interval& x) {
  x= Interval (Down::template euler<C>(), Up::template euler<C>()); }
TMPL inline void set_catalan (Interval& x) {
  x= Interval (Down::template catalan<C>(), Up::template catalan<C>()); }
TMPL inline Interval times_infinity (const Interval& x) {
  return x * infinity_cst<Interval > (); }

/******************************************************************************
* Basic arithmetic
******************************************************************************/

TMPL inline Interval
operator - (const Interval& x) {
  return Interval (-upper (x), -lower (x));
}

TMPL inline Interval
operator + (const Interval& x1, const Interval& x2) {
  return Interval (Down::add (lower (x1), lower (x2)),
                   Up::add (upper (x1), upper (x2)));
}

TMPL inline Interval
operator + (const Interval& x1, const C& x2) {
  return Interval (Down::add (lower (x1), x2), Up::add (upper (x1), x2));
}

TMPL inline Interval
operator + (const C& x1, const Interval& x2) {
  return Interval (Down::add (x1, lower (x2)), Up::add (x1, upper (x2)));
}

TMPL inline Interval
operator - (const Interval& x1, const Interval& x2) {
  return Interval (Down::sub (lower (x1), upper (x2)),
                   Up::sub (upper (x1), lower (x2)));
}

TMPL inline Interval
operator - (const Interval& x1, const C& x2) {
  return Interval (Down::add (lower (x1), x2), Up::add (upper (x1), x2));
}

TMPL inline Interval
operator - (const C& x1, const Interval& x2) {
  return Interval (Down::add (x1, upper (x2)), Up::add (x1, lower (x2)));
}

TMPL Interval
operator * (const Interval& x1, const Interval& x2) {
  if (lower (x1) >= 0)
    return Interval (lower (x2) >= 0?
                     Down::mul (lower (x1), lower (x2)):
                     Down::mul (upper (x1), lower (x2)),
                     upper (x2) >= 0?
                     Up::mul (upper (x1), upper (x2)):
                     Up::mul (lower (x1), upper (x2)));
  else if (upper (x1) <= 0)
    return Interval (upper (x2) >= 0?
                     Down::mul (lower (x1), upper (x2)):
                     Down::mul (upper (x1), upper (x2)),
                     lower (x2) >= 0?
                     Up::mul (upper (x1), lower (x2)):
                     Up::mul (lower (x1), lower (x2)));
  else if (lower (x2) >= 0)
    return Interval (Down::mul (lower (x1), upper (x2)),
                     Up::mul (upper (x1), upper (x2)));
  else if (upper (x2) <= 0)
    return Interval (Down::mul (upper (x1), lower (x2)),
                     Up::mul (lower (x1), lower (x2)));
  else
    return Interval (min (Down::mul (lower (x1), upper (x2)),
                          Down::mul (upper (x1), lower (x2))),
                     max (Up::mul (lower (x1), lower (x2)),
                          Up::mul (upper (x1), upper (x2))));
}

TMPL inline Interval operator * (const Interval& x1, const C& x2) {
  return x1 * Interval (x2); }
TMPL inline Interval operator * (const C& x1, const Interval& x2) {
  return Interval (x1) * x2; }

TMPL inline Interval
square (const Interval& x) {
  if (lower (x) >= 0)
    return Interval (Down::mul (lower (x), lower (x)),
                     Up::mul (upper (x), upper (x)));
  else if (upper (x) <= 0)
    return Interval (Down::mul (upper (x), upper (x)),
                     Up::mul (lower (x), lower (x)));
  else return Interval (Down::mul (lower (x), upper (x)),
                        max (Up::mul (lower (x), lower (x)),
                             Up::mul (upper (x), upper (x))));
}

TMPL Interval
operator / (const Interval& x1, const Interval& x2) {
  if (is_zero (x2)) return Nan (Interval);
  else if (lower (x2) > 0)
    return Interval (lower (x1) > 0?
                     Down::div (lower (x1), upper (x2)):
                     Down::div (lower (x1), lower (x2)),
                     upper (x1) > 0?
                     Up::div (upper (x1), lower (x2)):
                     Up::div (upper (x1), upper (x2)));
  else if (upper (x2) < 0)
    return Interval (upper (x1) > 0?
                     Down::div (upper (x1), upper (x2)):
                     Down::div (upper (x1), lower (x2)),
                     lower (x1) > 0?
                     Up::div (lower (x1), lower (x2)):
                     Up::div (lower (x1), upper (x2)));
  else return Nan (Interval);
}

TMPL inline Interval operator / (const Interval& x1, const C& x2) {
  return x1 / Interval (x2); }
TMPL inline Interval operator / (const C& x1, const Interval& x2) {
  return Interval (x1) / x2; }

ARITH_INT_SUGAR(TMPL,Interval);

/******************************************************************************
* In place arithmetic operations
******************************************************************************/

TMPL inline Interval&
Interval::operator += (const Interval& x) {
  l= Down::add (l, lower (x));
  r= Up::add (r, upper (x));
  return *this;
}

TMPL inline Interval&
Interval::operator -= (const Interval& x) {
  l= Down::sub (l, upper (x));
  r= Up::sub (r, lower (x));
  return *this;
}

/******************************************************************************
* Elementary functions
******************************************************************************/

TMPL Interval
sqrt (const Interval& x) {
  return Interval (Down::sqrt (lower (x)), Up::sqrt (upper (x)));
}

TMPL Interval
hypot (const Interval& x, const Interval& y) {
  Interval z= square (x) + square (y);
  if (z > 0) return sqrt (z);
  else return Interval (0, Up::sqrt (upper (z)));
}

TMPL Interval
exp (const Interval& x) {
  return Interval (Down::exp (lower (x)), Up::exp (upper (x)));
}

TMPL Interval
log (const Interval& x) {
  return Interval (Down::log (lower (x)), Up::log (upper (x)));
}

TMPL inline Interval
pow (const Interval& x, const Interval& y) {
  return exp (y * log (x));
}

TMPL Interval
cos (const Interval& x) {
  C pi_lo= promote (314159, lower (x)) / promote (100000, lower (x));
  if (Up::sub (upper (x), lower (x)) >= pi_lo)
    return Interval (promote (-1, lower (x)), promote (1, lower (x)));
  else {
    int s1= -sign (sin (lower (x)));
    int s2= -sign (sin (upper (x)));
    if (s1 >= 0 && s2 >= 0)
      return Interval (Down::cos (lower (x)), Up::cos (upper (x)));
    else if (s1 <= 0 && s2 <= 0)
      return Interval (Down::cos (upper (x)), Down::cos (lower (x)));
    else if (s1 >= 0)
      return Interval (min (Down::cos (lower (x)), Down::cos (upper (x))),
                       promote (1, lower (x)));
    else
      return Interval (promote (-1, lower (x)),
                       max (Down::cos (lower (x)), Down::cos (upper (x))));
  }
}

TMPL Interval
sin (const Interval& x) {
  C pi_lo= promote (314159, lower (x)) / promote (100000, lower (x));
  if (Up::sub (upper (x), lower (x)) >= pi_lo)
    return Interval (promote (-1, lower (x)), promote (1, lower (x)));
  else {
    int s1= sign (cos (lower (x)));
    int s2= sign (cos (upper (x)));
    if (s1 >= 0 && s2 >= 0)
      return Interval (Down::sin (lower (x)), Up::sin (upper (x)));
    else if (s1 <= 0 && s2 <= 0)
      return Interval (Down::sin (upper (x)), Down::sin (lower (x)));
    else if (s1 >= 0)
      return Interval (min (Down::sin (lower (x)), Down::sin (upper (x))),
                       promote (1, lower (x)));
    else
      return Interval (promote (-1, lower (x)),
                       max (Down::sin (lower (x)), Down::sin (upper (x))));
  }
}

TMPL Interval
tan (const Interval& x) {
  return sin (x) / cos (x);
}

TMPL inline Interval
cosh (const Interval& x) {
  if (is_nan (x)) return Nan (Interval);
  else if (lower (x) >= 0)
    return Interval (Down::cosh (lower (x)), Up::cosh (upper (x)));
  else if (upper (x) <= 0)
    return Interval (Down::cosh (upper (x)), Up::cosh (lower (x)));
  else return Interval (1, max (Up::cosh (lower (x)), Up::cosh (upper (x))));
}

TMPL inline Interval
sinh (const Interval& x) {
  return Interval (Down::sinh (lower (x)), Up::sinh (upper (x)));
}

TMPL Interval
tanh (const Interval& x) {
  return sinh (x) / cosh (x);
}

TMPL Interval
acos (const Interval& x) {
  return Interval (Down::acos (upper (x)), Up::acos (lower (x)));
}

TMPL Interval
asin (const Interval& x) {
  return Interval (Down::asin (lower (x)), Up::asin (upper (x)));
}

TMPL Interval
atan (const Interval& x) {
  return Interval (Down::atan (lower (x)), Up::atan (upper (x)));
}

TMPL Interval
atan2 (const Interval& y, const Interval& x) {
  if (sign (lower (x)) > 0) return atan (y/x);
  if (sign (lower (y)) > 0) return atan (-x/y) + (Pi (Interval) >> 1);
  if (sign (upper (y)) < 0) return atan (-x/y) - (Pi (Interval) >> 1);
  if (sign (upper (x)) < 0) return atan (y/x) + Pi (Interval);
  C pi= Up::template pi<C> ();
  return Interval (-pi, pi);
}

INV_TRIGO_SUGAR(TMPL,Interval)
INV_HYPER_SUGAR(TMPL,Interval)
ARG_HYPER_SUGAR(TMPL,Interval)

/******************************************************************************
* Floating point related functions
******************************************************************************/

TMPL inline bool is_finite (const Interval& x) {
  return is_finite (lower (x)) && is_finite (upper (x)); }
TMPL inline bool is_infinite (const Interval& x) {
  return is_infinite (lower (x)) && upper (x) == lower (x); }
TMPL inline bool is_fuzz (const Interval& x) {
  return is_infinite (lower (x)) && is_infinite (upper (x)) &&
    lower (x) != upper (x); }
TMPL inline bool is_nan (const Interval& x) {
  return is_nan (lower (x)) || is_nan (upper (x)); }
TMPL inline bool is_reliable (const Interval& x) {
  (void) x; return true; }

TMPL inline Interval change_precision (const Interval& x, xnat prec) {
  return Interval (change_precision (lower (x), prec),
                   change_precision (upper (x), prec)); }
TMPL inline xnat precision (const Interval& x) {
  return min (precision (lower (x)), precision (upper (x))); }

TMPL inline xint exponent (const Interval& x) {
  return max (exponent (lower (x)), exponent (upper (x))); }
TMPL inline double magnitude (const Interval& x) {
  return max (magnitude (lower (x)), magnitude (upper (x))); }
TMPL inline Interval operator << (const Interval& x, const xint& shift) {
  return Interval (incexp2 (lower (x), shift), incexp2 (upper (x), shift)); }
TMPL inline Interval operator >> (const Interval& x, const xint& shift) {
  return Interval (decexp2 (lower (x), shift), decexp2 (upper (x), shift)); }

TMPL inline Interval&
Interval::operator <<= (const xint& shift) {
  incexp2_assign (l, shift);
  incexp2_assign (r, shift);
  return *this;
}

TMPL inline Interval&
Interval::operator >>= (const xint& shift) {
  decexp2_assign (l, shift);
  decexp2_assign (r, shift);
  return *this;
}

TMPL inline Interval sharpen (const Interval& x) {
  return Interval (center (x)); }
TMPL inline Interval blur (const Interval& x, const C& r) {
  return x + Interval (-r, r); }
TMPL inline Interval blur (const Interval& x, const Interval& r) {
  return blur (x, abs_up (r)); }

#undef TMPL_DEF
#undef TMPL
#undef Interval
#undef Up
#undef Down
} // namespace mmx
#endif // __MMX_INTERVAL_HPP