# Interval Arithmetic Library

Contents of this page: Introduction Synopsis Template class `interval` Operations and functions Interval support library Compilation notes Common pitfalls and dangers Rationale History and Acknowledgments | Other pages associated with this page: Rounding policies Checking policies Policies manipulation Comparisons Base number type requirements Choosing your own interval type Test and example programs Headers inclusion Some items on the todo list

## Introduction and Overview

As implied by its name, this library is intended to help manipulating mathematical intervals. It consists of a single header <boost/numeric/interval.hpp> and principally a type which can be used as `interval<T>`. In fact, this interval template is declared as `interval<T,Policies>` where `Policies` is a policy class that controls the various behaviours of the interval class; `interval<T>` just happens to pick the default policies for the type `T`.

### Interval Arithmetic

An interval is a pair of numbers which represents all the numbers between these two. (Intervals are considered close so the bounds are included.) The purpose of this library is to extend the usual arithmetic functions to intervals. These intervals will be written [a,b] to represent all the numbers between a and b (included). a and b can be infinite (but they can not be the same infinite) and a ≤ b.

The fundamental property of interval arithmetic is the inclusion property:

  ``if f is a function on a set of numbers, f can be extended to a new function defined on intervals. This new function f takes one interval argument and returns an interval result such as: ∀ x ∈ [a,b], f(x) ∈ f([a,b]).''

Such a property is not limited to functions with only one argument. Whenever possible, the interval result should be the smallest one able to satisfy the property (it is not really useful if the new functions always answer [-∞,+∞]).

There are at least two reasons a user would like to use this library. The obvious one is when the user has to compute with intervals. One example is when input data have some builtin imprecision: instead of a number, an input variable can be passed as an interval. Another example application is to solve equations, by bisecting an interval until the interval width is small enough. A third example application is in computer graphics, where computations with boxes, segments or rays can be reduced to computations with points via intervals.

Another common reason to use interval arithmetic is when the computer doesn't produce exact results: by using intervals, it is possible to quantify the propagation of rounding errors. This approach is used often in numerical computation. For example, let's assume the computer stores numbers with ten decimal significant digits. To the question 1 + 1E-100 - 1, the computer will answer 0 although the correct answer would be 1E-100. With the help of interval arithmetic, the computer will answer [0,1E-9]. This is quite a huge interval for such a little result, but the precision is now known, without having to compute error propagation.

### Numbers, rounding, and exceptional behavior

The base number type is the type that holds the bounds of the interval. In order to successfully use interval arithmetic, the base number type must present some characteristics. Firstly, due to the definition of an interval, the base numbers have to be totally ordered so, for instance, `complex<T>` is not usable as base number type for intervals. The mathematical functions for the base number type should also be compatible with the total order (for instance if x>y and z>t, then it should also hold that x+z > y+t), so modulo types are not usable either.

Secondly, the computations must be exact or provide some rounding methods (for instance, toward minus or plus infinity) if we want to guarantee the inclusion property. Note that we also may explicitely specify no rounding, for instance if the base number type is exact, i.e. the result of a mathematic operation is always computed and represented without loss of precision. If the number type is not exact, we may still explicitely specify no rounding, with the obvious consequence that the inclusion property is no longuer guaranteed.

Finally, because heavy loss of precision is always possible, some numbers have to represent infinities or an exceptional behavior must be defined. The same situation also occurs for NaN (Not a Number).

Given all this, one may want to limit the template argument T of the class template `interval` to the floating point types `float`, `double`, and `long double`, as defined by the IEEE-754 Standard. Indeed, if the interval arithmetic is intended to replace the arithmetic provided by the floating point unit of a processor, these types are the best choice. Unlike `std::complex`, however, we don't want to limit `T` to these types. This is why we allow the rounding and exceptional behaviors to be given by the two policies (rounding and checking). We do nevertheless provide highly optimized rounding and checking class specializations for the above-mentioned floating point types.

### Operations and functions

It is straightforward to define the elementary arithmetic operations on intervals, being guided by the inclusion property. For instance, if [a,b] and [c,d] are intervals, [a,b]+[c,d] can be computed by taking the smallest interval that contains all the numbers x+y for x in [a,b] and y in [c,d]; in this case, rounding a+b down and c+d up will suffice. Other operators and functions are similarly defined (see their definitions below).

### Comparisons

It is also possible to define some comparison operators. Given two intervals, the result is a tri-state boolean type {false,true,indeterminate}. The answers false and true are easy to manipulate since they can directly be mapped on the boolean true and false. But it is not the case for the answer indeterminate since comparison operators are supposed to be boolean functions. So, what to do in order to obtain boolean answers?

One solution consists of deciding to adopt an exceptional behavior, such as a failed assertion or raising an exception. In this case, the exceptional behavior will be triggered when the result is indeterminate.

Another solution is to map indeterminate always to false, or always to true. If false is chosen, the comparison will be called "certain;" indeed, the result of [a,b] < [c,d] will be true if and only if: ∀ x ∈ [a,b] ∀ y ∈ [c,d], x < y. If true is chosen, the comparison will be called "possible;" indeed, the result of [a,b] < [c,d] will be true if and only if: ∃ x ∈ [a,b] ∃ y ∈ [c,d], x < y.

Since any of these solution has a clearly defined semantics, it is not clear that we should enforce either of them. For this reason, the default behavior consists to mimic the real comparisons by throwing an exception in the indeterminate case. Other behaviors can be selected bu using specific comparison namespace. There is also a bunch of explicitely named comparison functions. See comparisons pages for further details.

### Overview of the library, and usage

This library provides two quite distinct levels of usage. One is to use the basic class template `interval<T>` without specifying the policy. This only requires to know and understand the concepts developed above and the content of the namespace boost. In addition to the class `interval<T>`, this level of usage provides arithmetic operators (`+`, `-`, `*`, `/`), algebraic and piecewise-algebraic functions (`abs`, `square`, `sqrt`, `pow`), transcendental and trigonometric functions (`exp`, `log`, `sin`, `cos`, `tan`, `asin`, `acos`, `atan`, `sinh`, `cosh`, `tanh`, `asinh`, `acosh`, `atanh`), and the standard comparison operators (`<`, `<=`, `>`, `>=`, `==`, `!=`), as well as several interval-specific functions (`min`, `max`, which have a different meaning than `std::min` and `std::max`; `lower`, `upper`, `width`, `median`, `empty`, `singleton, equal`, `in`, `in_zero`, `subset`, `proper_subset`, `overlap`, `intersection`, `hull`, `bisect`).

For some functions which take several parameters of type `interval<T>`, all combinations of argument types `T` and `interval<T>` which contain at least one `interval<T>`, are considered in order to avoid a conversion from the arguments of type `T` to a singleton of type `interval<T>`. This is done for efficiency reasons (the fact that an argument is a singleton sometimes renders some tests unnecessary).

A somewhat more advanced usage of this library is to hand-pick the policies `Rounding` and `Checking` and pass them to `interval<T, Policies>` through the use of `Policies := boost::interval_lib::policies<Rounding,Checking>`. Appropriate policies can be fabricated by using the various classes provided in the namespace `boost::interval_lib` as detailed in section Interval Support Library. It is also possible to choose the comparison scheme by overloading operators through namespaces.

## Synopsis

    namespace boost {
    namespace numeric {

    namespace interval_lib {

    /* this declaration is necessary for the declaration of interval */
    template <class T> struct default_policies;

    /* ... ; the full synopsis of namespace interval_lib can be found here */

    } // namespace interval_lib

    /* template interval_policies; class definition can be found here */
    template<class Rounding, class Checking>
    struct interval_policies;

    /* template class interval; class definition can be found here */
    template<class T, class Policies = typename interval_lib::default_policies<T> >
    class interval;

    /* arithmetic operators involving intervals */
    template <class T, class Traits>  interval<T, Traits> operator+(const interval<T, Traits>& x);
    template <class T, class Traits>  interval<T, Traits> operator-(const interval<T, Traits>& x);

    template <class T, class Traits>  interval<T, Traits> operator+(const interval<T, Traits>& x, const interval<T, Traits>& y);
    template <class T, class Traits>  interval<T, Traits> operator+(const interval<T, Traits>& x, const T& y);
    template <class T, class Traits>  interval<T, Traits> operator+(const T& x, const interval<T, Traits>& y);

    template <class T, class Traits>  interval<T, Traits> operator-(const interval<T, Traits>& x, const interval<T, Traits>& y);
    template <class T, class Traits>  interval<T, Traits> operator-(const interval<T, Traits>& x, const T& y);
    template <class T, class Traits>  interval<T, Traits> operator-(const T& x, const interval<T, Traits>& y);

    template <class T, class Traits>  interval<T, Traits> operator*(const interval<T, Traits>& x, const interval<T, Traits>& y);
    template <class T, class Traits>  interval<T, Traits> operator*(const interval<T, Traits>& x, const T& y);
    template <class T, class Traits>  interval<T, Traits> operator*(const T& x, const interval<T, Traits>& y);

    template <class T, class Traits>  interval<T, Traits> operator/(const interval<T, Traits>& x, const interval<T, Traits>& y);
    template <class T, class Traits>  interval<T, Traits> operator/(const interval<T, Traits>& x, const T& y);
    template <class T, class Traits>  interval<T, Traits> operator/(const T& r, const interval<T, Traits>& x);

    /* algebraic functions: sqrt, abs, square, pow */
    template <class T, class Traits>  interval<T, Traits> abs(const interval<T, Traits>& x);
    template <class T, class Traits>  interval<T, Traits> sqrt(const interval<T, Traits>& x);
    template <class T, class Traits>  interval<T, Traits> square(const interval<T, Traits>& x);
    template <class T, class Traits>  interval<T, Traits> pow(const interval<T, Traits>& x, int y);

    /* transcendental functions: exp, log */
    template <class T, class Traits>  interval<T, Traits> exp(const interval<T, Traits>& x);
    template <class T, class Traits>  interval<T, Traits> log(const interval<T, Traits>& x);

    /* fmod, for trigonometric function argument reduction (see below) */
    template <class T, class Traits>  interval<T, Traits> fmod(const interval<T, Traits>& x, const interval<T, Traits>& y);
    template <class T, class Traits>  interval<T, Traits> fmod(const interval<T, Traits>& x, const T& y);
    template <class T, class Traits>  interval<T, Traits> fmod(const T& x, const interval<T, Traits>& y);

    /* trigonometric functions */
    template <class T, class Traits>  interval<T, Traits> sin(const interval<T, Traits>& x);
    template <class T, class Traits>  interval<T, Traits> cos(const interval<T, Traits>& x);
    template <class T, class Traits>  interval<T, Traits> tan(const interval<T, Traits>& x);
    template <class T, class Traits>  interval<T, Traits> asin(const interval<T, Traits>& x);
    template <class T, class Traits>  interval<T, Traits> acos(const interval<T, Traits>& x);
    template <class T, class Traits>  interval<T, Traits> atan(const interval<T, Traits>& x);

    /* hyperbolic trigonometric functions */
    template <class T, class Traits>  interval<T, Traits> sinh(const interval<T, Traits>& x);
    template <class T, class Traits>  interval<T, Traits> cosh(const interval<T, Traits>& x);
    template <class T, class Traits>  interval<T, Traits> tanh(const interval<T, Traits>& x);
    template <class T, class Traits>  interval<T, Traits> asinh(const interval<T, Traits>& x);
    template <class T, class Traits>  interval<T, Traits> acosh(const interval<T, Traits>& x);
    template <class T, class Traits>  interval<T, Traits> atanh(const interval<T, Traits>& x);

    /* min, max external functions (NOT std::min/max, see below) */
    template <class T, class Traits>  interval<T, Traits> max(const interval<T, Traits>& x, const interval<T, Traits>& y);
    template <class T, class Traits>  interval<T, Traits> max(const interval<T, Traits>& x, const T& y);
    template <class T, class Traits>  interval<T, Traits> max(const T& x, const interval<T, Traits>& y);
    template <class T, class Traits>  interval<T, Traits> min(const interval<T, Traits>& x, const interval<T, Traits>& y);
    template <class T, class Traits>  interval<T, Traits> min(const interval<T, Traits>& x, const T& y);
    template <class T, class Traits>  interval<T, Traits> min(const T& x, const interval<T, Traits>& y);

    /* bounds-related interval functions */
    template <class T, class Traits>  T lower(const interval<T, Traits>& x);
    template <class T, class Traits>  T upper(const interval<T, Traits>& x);
    template <class T, class Traits>  T width(const interval<T, Traits>& x);
    template <class T, class Traits>  T median(const interval<T, Traits>& x);

    /* bounds-related interval functions */
    template <class T, class Traits>  bool empty(const interval<T, Traits>& b);
    template <class T, class Traits>  bool singleton(const interval<T, Traits>& x);
    template <class T, class Traits>  bool equal(const interval<T, Traits>& x, const interval<T, Traits>& y);
    template <class T, class Traits>  bool in(const T& r, const interval<T, Traits>& b);
    template <class T, class Traits>  bool in_zero(const interval<T, Traits>& b);
    template <class T, class Traits>  bool subset(const interval<T, Traits>& a, const interval<T, Traits>& b);
    template <class T, class Traits>  bool proper_subset(const interval<T, Traits>& a, const interval<T, Traits>& b);
    template <class T, class Traits>  bool overlap(const interval<T, Traits>& x, const interval<T, Traits>& y);

    /* set manipulation interval functions */
    template <class T, class Traits>  interval<T, Traits> intersection(const interval<T, Traits>& x, const interval<T, Traits>& y);
    template <class T, class Traits>  interval<T, Traits> hull(const interval<T, Traits>& x, const interval<T, Traits>& y);
    template <class T, class Traits>  interval<T, Traits> hull(const interval<T, Traits>& x, const T& y);
    template <class T, class Traits>  interval<T, Traits> hull(const T& x, const interval<T, Traits>& y);
    template <class T, class Traits>  interval<T, Traits> hull(const T& x, const T& y);
    template <class T, class Traits>  std::pair<interval<T, Traits>, interval<T, Traits> > bisect(const interval<T, Traits>& x);

    /* interval comparison operators */
    template<class T, class Traits>  bool operator<(const interval<T, Traits>& x, const interval<T, Traits>& y);
    template<class T, class Traits>  bool operator<(const interval<T, Traits>& x, const T& y);
    template<class T, class Traits>  bool operator<(const T& x, const interval<T, Traits>& y);

    template<class T, class Traits>  bool operator<=(const interval<T, Traits>& x, const interval<T, Traits>& y);
    template<class T, class Traits>  bool operator<=(const interval<T, Traits>& x, const T& y);
    template<class T, class Traits>  bool operator<=(const T& x, const interval<T, Traits>& y);

    template<class T, class Traits>  bool operator>(const interval<T, Traits>& x, const interval<T, Traits>& y);
    template<class T, class Traits>  bool operator>(const interval<T, Traits>& x, const T& y);
    template<class T, class Traits>  bool operator>(const T& x, const interval<T, Traits>& y);

    template<class T, class Traits>  bool operator>=(const interval<T, Traits>& x, const interval<T, Traits>& y);
    template<class T, class Traits>  bool operator>=(const interval<T, Traits>& x, const T& y);
    template<class T, class Traits>  bool operator>=(const T& x, const interval<T, Traits>& y);

    template<class T, class Traits>  bool operator==(const interval<T, Traits>& x, const interval<T, Traits>& y);
    template<class T, class Traits>  bool operator==(const interval<T, Traits>& x, const T& y);
    template<class T, class Traits>  bool operator==(const T& x, const interval<T, Traits>& y);

    template<class T, class Traits>  bool operator!=(const interval<T, Traits>& x, const interval<T, Traits>& y);
    template<class T, class Traits>  bool operator!=(const interval<T, Traits>& x, const T& y);
    template<class T, class Traits>  bool operator!=(const T& x, const interval<T, Traits>& y);

    } // namespace numeric
    } // namespace boost

## Template class `interval`

The public interface of the template class interval itself is kept at a simplest minimum:

    template <class T, class Policies = typename interval_lib::default_policies<T>::type>
    class interval
    {
      public:
      typedef T      base_type;
      typedef Traits traits_type;

      interval(const T& v = T());
      interval(const T& l, const T& u);

      template<class Traits2>
      interval(const interval<T,Traits2>& r);

      // compiler-generated copy constructor and assignment operator are fine

      interval& operator=(const T& x);
      void assign(const T& l, const T& u);

      const T& lower() const;
      const T& upper() const;

      static interval empty();
      static interval whole();
      static interval hull(const T& x, const T& y);

      interval& operator+= (const T& r);
      interval& operator+= (const interval& r);
      interval& operator-= (const T& r);
      interval& operator-= (const interval& r);
      interval& operator*= (const T& r);
      interval& operator*= (const interval& r);
      interval& operator/= (const T& r);
      interval& operator/= (const interval& r);
    };

The two main constructors accept one and two arguments of type T. The first one produces a singleton interval [v,v]; and the second one an interval [l,u]. The second one will use the checking policy if u<l. If you are not sure whether the bounds are inverted or not, it is better to use the function call `hull(l,u)`. This will not create a problem if u<l.

There is a template constructor in order to change the traits parameter of an interval. But there is no constructor that changes the base type of an interval (it could lead to violating the inclusion property if the old base type is not contained in the new one).

There is an assign function in order to directly change the bounds of an interval. It behaves like the second constructor if the bounds are inverted. There is no assign function that directly takes an interval or only one number as a parameter; just use the assignment operator in that case.

Assignment operator, copy constructor and destructor are the ones the compiler automatically generates. There is also an assignment operator for the type T.

The static functions `empty` and `whole` produces the corresponding intervals. They are static member functions rather than global functions because they cannot guess their return types. Likewise for `hull`.

## Operations and Functions

Some of the following functions expect `min` and `max` to be defined for the base type. Those are the only requirements for the `interval` class (but the policies can have other requirements).

#### Operators `+` `-` `*` `/` `+=` `-=` `*=` `/=`

The basic operations are the unary minus and the binary `+` `-` `*` `/`. The unary minus takes an interval and returns an interval. The binary operations take two intervals, or one interval and a number, and return an interval. If an argument is a number instead of an interval, you can expect the result to be the same as if the number was first converted to an interval. This property will be true for all the following functions and operators.

There are also some assignment operators `+=` `-=` `*=` `/=`. There is not much to say: `x op= y` is equivalent to `x = x op y`. If an exception is thrown during the computations, the l-value is not modified (but it may be corrupt if an exception is thrown by the base type during an assignment).

The operators `/` and `/=` will try to produce an empty interval if the denominator is exactly zero. If the denominator contains zero (but not only zero), the result will be the smallest interval containing the set of division results; so one of its bound will be infinite, but it may not be the whole interval.

#### `lower` `upper` `median`

`lower`, `upper`, `median` respectively compute the lower bound, the upper bound, and the median number of an interval (`(lower+upper)/2` rounded to nearest). `width` computes the width of an interval (`upper-lower` rounded toward plus infinity).

#### `min` `max` `abs` `square` `pow` `division_part?`

The functions `min`, `max` and `abs` are also defined. Please do not mistake them for the functions defined in the standard library (aka `a<b?a:b`, `a>b?a:b`, `a<0?-a:a`). These functions are compatible with the elementary property of interval arithmetic. For example, max([a,b], [c,d]) = {max(x,y) such that x in [a,b] and y in [c,d]}. They are not defined in the `std` namespace but in the boost namespace in order to avoid conflict with the other definitions.

The `square` function is quite particular. As you can expect from its name, it computes the square of its argument. The reason this function is provided is: `square(x)` is not `x*x` but only a subset when `x` contains zero. For example, [-2,2]*[-2,2] = [-4,4] but [-2,2]˛ = [0,4]; the result is a smaller interval. Consequently, `square(x)` should be used instead of `x*x` because of its better accuracy and a small performance improvement.

As for `square`, `pow` provides an efficient and more accurate way to compute the integer power of an interval. Please note: when the power is 0 and the interval is not empty, the result is 1, even if the input interval contains 0.

The functions `division_part1` and `division_part2` are useful when the user expects the division to return disjoint intervals if necessary. For example, the narrowest closed set containg [2,3] / [-2,1] is not ]-∞,∞[ but the union of ]-∞,-1] and [2,∞[. When the result of the division is representable by only one interval, `division_part1` returns this interval and sets the boolean reference to `false`. However, if the result needs two intervals, `division_part1` returns the negative part and sets the boolean reference to `true`; a call to `division_part2` is now needed to get the positive part.

#### `intersect` `hull` `overlap` `in` `in_zero` `subset` `proper_subset` `empty` `singleton` `equal`

`intersect` computes the set intersection of two closed sets, `hull` computes the smallest interval which contains the two parameters; those parameters can be numbers or intervals. If one of the arguments is an invalid number or an empty interval, the function will only use the other argument to compute the resulting interval (if allowed by the checking policy).

There is no union function since the union of two intervals is not an interval if they do not overlap. If they overlap, the `hull` function computes the union.

The function `overlap` tests if two intervals have some common subset. `in` tests if a number is in an interval; `in_zero` is a variant which tests if zero is in the interval. `subset` tests if the first interval is a subset of the second; and `proper_subset` tests if it is a proper subset. `empty` and `singleton` test if an interval is empty or is a singleton. Finally, `equal` tests if two intervals are equal.

#### `sqrt` `log` `exp` `sin` `cos` `tan` `asin` `acos` `atan` `sinh` `cosh` `tanh` `asinh` `acosh` `atanh` `fmod`

The functions `sqrt`, `log`, `exp`, `sin`, `cos`, `tan`, `asin`, `acos`, `atan`, `sinh`, `cosh`, `tanh`, `asinh`, `acosh`, `atanh` are also defined. There is not much to say; these functions extend the traditional functions to the intervals and respect the basic property of interval arithmetic. They use the checking policy to produce empty intervals when the input interval is strictly outside of the domain of the function.

The function `fmod(interval x, interval y)` expects the lower bound of `y` to be strictly positive and returns an interval `z` such as `0 <= z.lower() < y.upper()` and such as `z` is a superset of `x-n*y` (with `n` being an integer). So, if the two arguments are positive singletons, this function `fmod(interval, interval)` will behave like the traditional function `fmod(double, double)`.

Please note that `fmod` does not respect the inclusion property of arithmetic interval. For example, the result of `fmod`([13,17],[7,8]) should be [0,8] (since it must contain [0,3] and [5,8]). But this answer is not really useful when the purpose is to restrict an interval in order to compute a periodic function. It is the reason why `fmod` will answer [5,10].

### Constants

Some constants are hidden in the `boost::interval_lib` namespace. They need to be explicitely templated by the interval type. The functions are `pi<I>()`, `pi_half<I>()` and `pi_twice<I>()`, and they return an object of interval type `I`. Their respective values are π, π/2 and 2π.

### Exception throwing

The interval class and all the functions defined around this class never throw any exceptions by themselves. However, it does not mean that an operation will never throw an exception. For example, let's consider the copy constructor. As explained before, it is the default copy constructor generated by the compiler. So it will not throw an exception if the copy constructor of the base type does not throw an exception.

The same situation applies to all the functions: exceptions will only be thrown if the base type or one of the three policies throws an exception.

## Interval Support Library

The interval support library consists of a collection of classes that can be used and combined to fabricate almost various commonly-needed interval policies. In contrast to the basic classes and functions which are used in conjunction with `interval<T>` (and the default policies as the implicit second template parameter in this type), which belong simply to the namespace `boost`, these components belong to the namespace `boost::numeric::interval_lib`.

We merely give the synopsis here and defer each section to a separate web page since it is only intended for the advanced user. This allows to expand on each topic with examples, without unduly stretching the limits of this document.

#### Synopsis

    namespace boost {
    namespace numeric {
    namespace interval_lib {

    /* built-in rounding policy and its specializations */
    template <class T>  struct rounded_math;
    template <>         struct rounded_math<float>;
    template <>         struct rounded_math<double>;
    template <>         struct rounded_math<long double>;

    /* built-in rounding construction blocks */
    template <class T>  struct rounding_control;

    template <class T, class Rounding = rounding_control<T> >  struct rounded_arith_exact;
    template <class T, class Rounding = rounding_control<T> >  struct rounded_arith_std;
    template <class T, class Rounding = rounding_control<T> >  struct rounded_arith_opp;

    template <class T, class Rounding>  struct rounded_transc_dummy;
    template <class T, class Rounding = rounded_arith_exact<T> >  struct rounded_transc_exact;
    template <class T, class Rounding = rounded_arith_std<T> >  struct rounded_transc_std;
    template <class T, class Rounding = rounded_arith_opp<T> >  struct rounded_transc_opp;

    template <class Rounding> struct save_state;
    template <class Rounding> struct save_state_nothing;

    /* built-in checking policies */
    template <class T> struct checking_nothing;
    template <class T> struct checking_lax;
    template <class T> struct checking_strict;

    /* some metaprogramming to manipulate interval policies */
    template <class Rounding, class Checking> struct policies;
    template <class OldInterval, class NewRounding> struct change_rounding;
    template <class OldInterval, class NewChecking> struct change_checking;
    template <class OldInterval> struct unprotect;

    /* constants, need to be explicitly templated */
    template<class I> I pi();
    template<class I> I pi_half();
    template<class I> I pi_twice();

    /* interval explicit comparison functions:
     * the mode can be cer=certainly or pos=possibly,
     * the function lt=less_than, gt=greater_than, le=less_than_or_equal_to, ge=greater_than_or_equal_to
     *   eq=equal_to, ne= not_equal_to */
    template <class T, class Traits>  bool cerlt(const interval<T, Traits>& x, const interval<T, Traits>& y);
    template <class T, class Traits>  bool cerlt(const interval<T, Traits>& x, const T& y);
    template <class T, class Traits>  bool cerlt(const T& x, const interval<T, Traits>& y);

    template <class T, class Traits>  bool cerle(const interval<T, Traits>& x, const interval<T, Traits>& y);
    template <class T, class Traits>  bool cerle(const interval<T, Traits>& x, const T& y);
    template <class T, class Traits>  bool cerle(const T& x, const interval<T, Traits>& y);

    template <class T, class Traits>  bool cergt(const interval<T, Traits>& x, const interval<T, Traits>& y);
    template <class T, class Traits>  bool cergt(const interval<T, Traits>& x, const T& y);
    template <class T, class Traits>  bool cergt(const T& x, const interval<T, Traits>& y);

    template <class T, class Traits>  bool cerge(const interval<T, Traits>& x, const interval<T, Traits>& y);
    template <class T, class Traits>  bool cerge(const interval<T, Traits>& x, const T& y);
    template <class T, class Traits>  bool cerge(const T& x, const interval<T, Traits>& y);

    template <class T, class Traits>  bool cereq(const interval<T, Traits>& x, const interval<T, Traits>& y);
    template <class T, class Traits>  bool cereq(const interval<T, Traits>& x, const T& y);
    template <class T, class Traits>  bool cereq(const T& x, const interval<T, Traits>& y);

    template <class T, class Traits>  bool cerne(const interval<T, Traits>& x, const interval<T, Traits>& y);
    template <class T, class Traits>  bool cerne(const interval<T, Traits>& x, const T& y);
    template <class T, class Traits>  bool cerne(const T& x, const interval<T, Traits>& y);

    template <class T, class Traits>  bool poslt(const interval<T, Traits>& x, const interval<T, Traits>& y);
    template <class T, class Traits>  bool poslt(const interval<T, Traits>& x, const T& y);
    template <class T, class Traits>  bool poslt(const T& x, const interval<T, Traits>& y);

    template <class T, class Traits>  bool posle(const interval<T, Traits>& x, const interval<T, Traits>& y);
    template <class T, class Traits>  bool posle(const interval<T, Traits>& x, const T& y);
    template <class T, class Traits>  bool posle(const T& x, const interval<T, Traits>& y);

    template <class T, class Traits>  bool posgt(const interval<T, Traits>& x, const interval<T, Traits>& y);
    template <class T, class Traits>  bool posgt(const interval<T, Traits>& x, const T& y);
    template <class T, class Traits>  bool posgt(const T& x, const interval<T, Traits> & y);

    template <class T, class Traits>  bool posge(const interval<T, Traits>& x, const interval<T, Traits>& y);
    template <class T, class Traits>  bool posge(const interval<T, Traits>& x, const T& y);
    template <class T, class Traits>  bool posge(const T& x, const interval<T, Traits>& y);

    template <class T, class Traits>  bool poseq(const interval<T, Traits>& x, const interval<T, Traits>& y);
    template <class T, class Traits>  bool poseq(const interval<T, Traits>& x, const T& y);
    template <class T, class Traits>  bool poseq(const T& x, const interval<T, Traits>& y);

    template <class T, class Traits>  bool posne(const interval<T, Traits>& x, const interval<T, Traits>& y);
    template <class T, class Traits>  bool posne(const interval<T, Traits>& x, const T& y);
    template <class T, class Traits>  bool posne(const T& x, const interval<T, Traits>& y);

    /* comparison namespaces */
    namespace compare {
      namespace lexicographic;
      namespace set;
      namespace tribool;
    } // namespace compare

    } // namespace interval_lib
    } // namespace numeric
    } // namespace boost

Each component of the interval support library is detailed in its own page.

- Comparisons
- Rounding
- Checking

## Compilation Notes

#### Compilation options

We have made many efforts to reduce the number of compilation options to a minimum, but depending on your platform, you may find that you need one or two. Some versions of gcc miss proper specialization of `numeric_limits` and using the flag `-DBOOST_NO_LIMITS` may help. Also removing the warnings for older gcc compilers about the lack of long double support can be done with the option `-Wno-long-double`. The tests and examples have successfully been compiled with gcc 2.95 on x86 (Debian, FreeBSD), on sparc (Debian); with gcc 2.95 and stlport on x86 (Solaris), sparc (Solaris), and powerpc (Debian); with gcc 2.96, 3.04, 3.1.1, and 3.2 on x86 (Debian); with Intel cc 6.0 and 7.0 on x86 (Debian); and with SunPRO 5.3 (Solaris). There are some issues with Borland on x86 (Windows), and Visual C++ which we're working on.

## Common Pitfalls and Dangers

#### Comparisons

One of the biggest problems is problably the correct use of the comparison functions and operators. First, functions and operators do not try to know if two intervals are the same mathematical object. So, if the comparison used is "certain", then `x != x` is always true unless `x` is a singleton interval; and the same problem arises with `cereq` and `cerne`.

Another misleading interpretation of the comparison is: you cannot always expect [a,b] < [c,d] to be !([a,b] >= [c,d]) since the comparison is not necessarily total. Equality and less comparison should be seen as two distincts relational operators.

#### Interval values and references

This problem is a corollary of the previous problem with `x != x`. All the functions of the library only consider the value of an interval and not the reference of an interval. In particular, you should not expect (unless `x` is a singleton) the following values to be equal: `x/x` and 1, `x*x` and `square(x)`, `x-x` and 0, etc. So the main cause of wide intervals is that interval arithmetic does not identify different occurences of the same variable. So, whenever possible, the user has to rewrite the formulas to eliminate multiple occurences of the same variable. For example, `square(x)-2*x` is far less precise than `square(x-1)-1`.

#### Unprotected rounding

As explained in this section, a good way to speed up computations when the base type is a basic floating-point type is to unprotect the intervals at the hot spots of the algorithm. This method is safe and really an improvement for interval computations. But please remember that any basic floating-point operation executed inside the unprotection blocks will probably have an undefined behavior (but only for the current thread). And do not forget to create a rounding object as explained in the example.

## Rationale

The purpose of this library is to provide an efficient and generalized way to deal with interval arithmetic through the use of a templatized class `boost::interval`. The big contention for which we provide a rationale is the format of this class template.

It would have been easier to provide a class interval whose base type is double. Or to follow `std::complex` and allow only specializations for `float`, `double`, and `long double`. We decided not to do this to allow intervals on custom types, e.g. fixed-precision bigfloat library types (MPFR, etc.)

Policy design. Although it was tempting to make it a class template with only one template argument, the diversity of uses for an interval arithmetic practically forced us to use policies. The behavior of this class can be fixed by three policies. These policies are packaged into a single policy class, rather than making `interval` with three template parameters. This is both for ease of use (the policy class can be picked by default) and for readability.

The first policy provides all the mathematical functions on the base type needed to define the functions on the interval type. The second one sets the way exceptional cases encountered during computations are handled.

We could foresee situations where any combination of these policies would be appropriate. Moreover, we wanted to enable the user of the library to reuse the `interval` class template while at the same time choosing his own behavior. See this page for some examples.

Rounding policy. The library provides specialized implementations of the second policy for the primitive types float and double. In order for these implementations to be correct and fast, the library needs to work a lot with rounding modes. Some processors are directly dealt with and some mecanisms are provided in order to speed up the computations. It seems to be heavy and hazardous optimizations for a gain of only a few computer cycles; but in reality, the speed-up factor can easily go past 2 or 3 depending on the computer. Moreover, these optimizations do not impact the interface in any major way (with the design we have chosen, everything can be added by specialization or by passing different template parameters).

Pred/succ. In a previous version, two functions `pred` and `succ`, with various corollaries like `widen`, were supplied. The intent was to enlarge the interval by one ulp (as little as possible), e.g. to ensure the inclusion property. Since making interval a template of T, we could not define ulp for a random parameter. In turn, rounding policies let us eliminate entirely the use of ulp while making the intervals tighter (if a result is a representable singleton, there is no use to widen the interval). We decided to drop those functions.

Specialization of `std::less`. Since the operator `<` depends on the comparison namespace chosen locally chosen by the user, it is not possible to correctly specialize `std::less`. So you have to explicitely provide such a class to all the algorithms and templates that could require it (for example, `std::map`).

Input/output. The basic header `boost/numeric/interval.hpp` does not include the I/O operators for the class interval. Somewhat default operators can be had by including `boost/numeric/interval/io.hpp`, and they make very few effort to guarantee the inclusion property (they do so for `float`, `double`, and `long double`). But they cannot guarantee it for any base number type, especially since it is not know how the type will be displayed (for example, there could be trouble if the type is a symbolic type that uses the character `[` or `]'). For this, we would need an extra i/o policy. It was not felt worth the trouble. For one thing, it's always possible to display the bounds. For another, intervals are used somewhat internally and there is seldom need to display them. But if this should be needed, the header is there.

Conversion from `interval<T1>` to `interval<T>`. There is no such conversion in the library. Specifically, there is no way to change the base number type without explicitly casting both the lower and upper bounds from `T1` to `T2` and using a constructor. We did not want to provide this conversion because there is no easy way to make sure the inclusion property is guaranteed: a simple casting with loss of precision can lead to invalid or empty intervals. Moreover it is not clear who should be responsible for the conversion: `T2` is not required to know about all the possible types `T1`. We also do not want to throw an exception is there is a loss of precision (as is done for instance in `numeric_cast`); instead, the result should be correctly rounded to ensure the inclusion property. So we decided not to provide the conversion and leave it up to the user to do it himself or herself.

## History and Acknowledgments

This library was mostly inspired by previous work from Jens Maurer. Some discussions about his work are reproduced here and the work itself can be found here. Jeremy Siek and Maarten Keijzer provided some rounding control for MSVC and Sparc platforms.

Guillaume Melquiond, Hervé Brönnimann and Sylvain Pion started from the library left by Jens and added the policy design. Guillaume and Sylvain worked hard on the code, especially the porting and mostly tuning of the rounding modes to the different architectures. Guillaume mostly did most of the coding, while Sylvain and Hervé have provided some useful comments in order for this library to be written. Hervé reorganized and wrote chapters of the documentation based on Guillaume's great starting point.

Revised: 2003-02-05
Copyright (c) Guillaume Melquiond, Sylvain Pion, Hervé Brönnimann, 2002-2003.
Polytechnic University, 2002.