Interval Arithmetic Library
Interval Arithmetic LibraryAs 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.
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:
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.
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.
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+c down and b+d up will suffice. Other operators and functions are similarly defined (see their definitions below).
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.
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::numeric::interval_lib::policies<Rounding,Checking>.
Appropriate policies can be fabricated by using the various classes provided
in the namespace boost::numeric::interval_lib as detailed in
section Interval Support Library. It is also
possible to choose the comparison scheme by overloading operators through
namespaces.
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>::type > class interval;
/* arithmetic operators involving intervals */
template <class T, class Policies> interval<T, Policies> operator+(const interval<T, Policies>& x);
template <class T, class Policies> interval<T, Policies> operator-(const interval<T, Policies>& x);
template <class T, class Policies> interval<T, Policies> operator+(const interval<T, Policies>& x, const interval<T, Policies>& y);
template <class T, class Policies> interval<T, Policies> operator+(const interval<T, Policies>& x, const T& y);
template <class T, class Policies> interval<T, Policies> operator+(const T& x, const interval<T, Policies>& y);
template <class T, class Policies> interval<T, Policies> operator-(const interval<T, Policies>& x, const interval<T, Policies>& y);
template <class T, class Policies> interval<T, Policies> operator-(const interval<T, Policies>& x, const T& y);
template <class T, class Policies> interval<T, Policies> operator-(const T& x, const interval<T, Policies>& y);
template <class T, class Policies> interval<T, Policies> operator*(const interval<T, Policies>& x, const interval<T, Policies>& y);
template <class T, class Policies> interval<T, Policies> operator*(const interval<T, Policies>& x, const T& y);
template <class T, class Policies> interval<T, Policies> operator*(const T& x, const interval<T, Policies>& y);
template <class T, class Policies> interval<T, Policies> operator/(const interval<T, Policies>& x, const interval<T, Policies>& y);
template <class T, class Policies> interval<T, Policies> operator/(const interval<T, Policies>& x, const T& y);
template <class T, class Policies> interval<T, Policies> operator/(const T& r, const interval<T, Policies>& x);
/* algebraic functions: sqrt, abs, square, pow */
template <class T, class Policies> interval<T, Policies> abs(const interval<T, Policies>& x);
template <class T, class Policies> interval<T, Policies> sqrt(const interval<T, Policies>& x);
template <class T, class Policies> interval<T, Policies> square(const interval<T, Policies>& x);
template <class T, class Policies> interval<T, Policies> pow(const interval<T, Policies>& x, int y);
/* transcendental functions: exp, log */
template <class T, class Policies> interval<T, Policies> exp(const interval<T, Policies>& x);
template <class T, class Policies> interval<T, Policies> log(const interval<T, Policies>& x);
/* fmod, for trigonometric function argument reduction (see below) */
template <class T, class Policies> interval<T, Policies> fmod(const interval<T, Policies>& x, const interval<T, Policies>& y);
template <class T, class Policies> interval<T, Policies> fmod(const interval<T, Policies>& x, const T& y);
template <class T, class Policies> interval<T, Policies> fmod(const T& x, const interval<T, Policies>& y);
/* trigonometric functions */
template <class T, class Policies> interval<T, Policies> sin(const interval<T, Policies>& x);
template <class T, class Policies> interval<T, Policies> cos(const interval<T, Policies>& x);
template <class T, class Policies> interval<T, Policies> tan(const interval<T, Policies>& x);
template <class T, class Policies> interval<T, Policies> asin(const interval<T, Policies>& x);
template <class T, class Policies> interval<T, Policies> acos(const interval<T, Policies>& x);
template <class T, class Policies> interval<T, Policies> atan(const interval<T, Policies>& x);
/* hyperbolic trigonometric functions */
template <class T, class Policies> interval<T, Policies> sinh(const interval<T, Policies>& x);
template <class T, class Policies> interval<T, Policies> cosh(const interval<T, Policies>& x);
template <class T, class Policies> interval<T, Policies> tanh(const interval<T, Policies>& x);
template <class T, class Policies> interval<T, Policies> asinh(const interval<T, Policies>& x);
template <class T, class Policies> interval<T, Policies> acosh(const interval<T, Policies>& x);
template <class T, class Policies> interval<T, Policies> atanh(const interval<T, Policies>& x);
/* min, max external functions (NOT std::min/max, see below) */
template <class T, class Policies> interval<T, Policies> max(const interval<T, Policies>& x, const interval<T, Policies>& y);
template <class T, class Policies> interval<T, Policies> max(const interval<T, Policies>& x, const T& y);
template <class T, class Policies> interval<T, Policies> max(const T& x, const interval<T, Policies>& y);
template <class T, class Policies> interval<T, Policies> min(const interval<T, Policies>& x, const interval<T, Policies>& y);
template <class T, class Policies> interval<T, Policies> min(const interval<T, Policies>& x, const T& y);
template <class T, class Policies> interval<T, Policies> min(const T& x, const interval<T, Policies>& y);
/* bounds-related interval functions */
template <class T, class Policies> T lower(const interval<T, Policies>& x);
template <class T, class Policies> T upper(const interval<T, Policies>& x);
template <class T, class Policies> T width(const interval<T, Policies>& x);
template <class T, class Policies> T median(const interval<T, Policies>& x);
template <class T, class Policies> T norm(const interval<T, Policies>& x);
/* bounds-related interval functions */
template <class T, class Policies> bool empty(const interval<T, Policies>& b);
template <class T, class Policies> bool singleton(const interval<T, Policies>& x);
template <class T, class Policies> bool equal(const interval<T, Policies>& x, const interval<T, Policies>& y);
template <class T, class Policies> bool in(const T& r, const interval<T, Policies>& b);
template <class T, class Policies> bool in_zero(const interval<T, Policies>& b);
template <class T, class Policies> bool subset(const interval<T, Policies>& a, const interval<T, Policies>& b);
template <class T, class Policies> bool proper_subset(const interval<T, Policies>& a, const interval<T, Policies>& b);
template <class T, class Policies> bool overlap(const interval<T, Policies>& x, const interval<T, Policies>& y);
/* set manipulation interval functions */
template <class T, class Policies> interval<T, Policies> intersection(const interval<T, Policies>& x, const interval<T, Policies>& y);
template <class T, class Policies> interval<T, Policies> hull(const interval<T, Policies>& x, const interval<T, Policies>& y);
template <class T, class Policies> interval<T, Policies> hull(const interval<T, Policies>& x, const T& y);
template <class T, class Policies> interval<T, Policies> hull(const T& x, const interval<T, Policies>& y);
template <class T, class Policies> interval<T, Policies> hull(const T& x, const T& y);
template <class T, class Policies> std::pair<interval<T, Policies>, interval<T, Policies> > bisect(const interval<T, Policies>& x);
/* interval comparison operators */
template<class T, class Policies> bool operator<(const interval<T, Policies>& x, const interval<T, Policies>& y);
template<class T, class Policies> bool operator<(const interval<T, Policies>& x, const T& y);
template<class T, class Policies> bool operator<(const T& x, const interval<T, Policies>& y);
template<class T, class Policies> bool operator<=(const interval<T, Policies>& x, const interval<T, Policies>& y);
template<class T, class Policies> bool operator<=(const interval<T, Policies>& x, const T& y);
template<class T, class Policies> bool operator<=(const T& x, const interval<T, Policies>& y);
template<class T, class Policies> bool operator>(const interval<T, Policies>& x, const interval<T, Policies>& y);
template<class T, class Policies> bool operator>(const interval<T, Policies>& x, const T& y);
template<class T, class Policies> bool operator>(const T& x, const interval<T, Policies>& y);
template<class T, class Policies> bool operator>=(const interval<T, Policies>& x, const interval<T, Policies>& y);
template<class T, class Policies> bool operator>=(const interval<T, Policies>& x, const T& y);
template<class T, class Policies> bool operator>=(const T& x, const interval<T, Policies>& y);
template<class T, class Policies> bool operator==(const interval<T, Policies>& x, const interval<T, Policies>& y);
template<class T, class Policies> bool operator==(const interval<T, Policies>& x, const T& y);
template<class T, class Policies> bool operator==(const T& x, const interval<T, Policies>& y);
template<class T, class Policies> bool operator!=(const interval<T, Policies>& x, const interval<T, Policies>& y);
template<class T, class Policies> bool operator!=(const interval<T, Policies>& x, const T& y);
template<class T, class Policies> bool operator!=(const T& x, const interval<T, Policies>& y);
namespace interval_lib {
template<class T, class Policies> interval<T, Policies> division_part1(const interval<T, Policies>& x, const interval<T, Policies& y, bool& b);
template<class T, class Policies> interval<T, Policies> division_part2(const interval<T, Policies>& x, const interval<T, Policies& y, bool b = true);
template<class T, class Policies> interval<T, Policies> multiplicative_inverse(const interval<T, Policies>& x);
template<class I> I add(const typename I::base_type& x, const typename I::base_type& y);
template<class I> I sub(const typename I::base_type& x, const typename I::base_type& y);
template<class I> I mul(const typename I::base_type& x, const typename I::base_type& y);
template<class I> I div(const typename I::base_type& x, const typename I::base_type& y);
} // namespace interval_lib
} // namespace numeric
} // namespace boost
intervaltemplate <class T, class Policies = typename interval_lib::default_policies<T>::type>
class interval
{
public:
typedef T base_type;
typedef Policies traits_type;
interval();
interval(T const &v);
template<class T1> interval(T1 const &v);
interval(T const &l, T const &u);
template<class T1, class T2> interval(T1 const &l, T2 const &u);
interval(interval<T, Policies> const &r);
template<class Policies1> interval(interval<T, Policies1> const &r);
template<class T1, class Policies1> interval(interval<T1, Policies1> const &r);
interval &operator=(T const &v);
template<class T1> interval &operator=(T1 const &v);
interval &operator=(interval<T, Policies> const &r);
template<class Policies1> interval &operator=(interval<T, Policies1> const &r);
template<class T1, class Policies1> interval &operator=(interval<T1, Policies1> const &r);
void assign(T const &l, T const &u);
T const &lower() const;
T const &upper() const;
static interval empty();
static interval whole();
static interval hull(T const &x, T const &y);
interval& operator+= (T const &r);
interval& operator-= (T const &r);
interval& operator*= (T const &r);
interval& operator/= (T const &r);
interval& operator+= (interval const &r);
interval& operator-= (interval const &r);
interval& operator*= (interval const &r);
interval& operator/= (interval const &r);
};
The constructors create an interval enclosing their arguments. If there are two arguments, the first one is assumed to be the left bound and the second one is the right bound. Consequently, the arguments need to be ordered. If the property !(l <= u) is not respected, the checking policy will be used to create an empty interval. If no argument is given, the created interval is the singleton zero.
If the type of the arguments is the same as the base number type, the
values are directly used for the bounds. If it is not the same type, the
library will use the rounding policy in order to convert the arguments
(conv_down and conv_up) and create an enclosing
interval. When the argument is an interval with a different policy, the input
interval is checked in order to correctly propagate its emptiness (if
empty).
The assignment operators behave similarly, except they obviously take one
argument only. There is also an assign function in order to
directly change the bounds of an interval. It behaves like the two-arguments
constructors if the bounds are not ordered. There is no assign function that
directly takes an interval or only one number as a parameter; just use the
assignment operators in such a case.
The static functions empty and whole produce the
corresponding intervals. They are static member functions rather than global
functions because they cannot guess their return types. Likewise for
hull. empty and whole involve the
checking policy in order to get the bounds of the resulting intervals.
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).
+ - * /
+= -= *= /=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
width normlower, 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). norm
computes an upper bound of the interval in absolute value; it is a
mathematical norm (hence the name) similar to the absolute value for real
numbers.
min max abs square
pow division_part?
multiplicative_inverseThe 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. multiplicative_inverse computes
1/x.
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. This second function can take the boolean
returned by the first function as last argument. If this bool is not given,
its value is assumed to be true and the behavior of the function is then
undetermined if the division does not produce a second interval.
intersect hull overlap
in in_zero subset
proper_subset empty singleton
equalintersect 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
fmodThe 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].
add sub mul divThese four functions take two numbers and return the enclosing interval for the operations. It avoids converting a number to an interval before an operation, it can result in a better code with poor optimizers.
Some constants are hidden in the boost::numeric::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π.
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 two policies throws an exception.
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.
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_base;
template <class T, class Checking = checking_base<T>, class Exception = exception_create_empty> struct checking_no_empty;
template <class T, class Checking = checking_base<T> > struct checking_no_nan;
template <class T, class Checking = checking_base<T>, class Exception = exception_invalid_number> struct checking_catch_nan;
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 Policies> bool cerlt(const interval<T, Policies>& x, const interval<T, Policies>& y);
template <class T, class Policies> bool cerlt(const interval<T, Policies>& x, const T& y);
template <class T, class Policies> bool cerlt(const T& x, const interval<T, Policies>& y);
template <class T, class Policies> bool cerle(const interval<T, Policies>& x, const interval<T, Policies>& y);
template <class T, class Policies> bool cerle(const interval<T, Policies>& x, const T& y);
template <class T, class Policies> bool cerle(const T& x, const interval<T, Policies>& y);
template <class T, class Policies> bool cergt(const interval<T, Policies>& x, const interval<T, Policies>& y);
template <class T, class Policies> bool cergt(const interval<T, Policies>& x, const T& y);
template <class T, class Policies> bool cergt(const T& x, const interval<T, Policies>& y);
template <class T, class Policies> bool cerge(const interval<T, Policies>& x, const interval<T, Policies>& y);
template <class T, class Policies> bool cerge(const interval<T, Policies>& x, const T& y);
template <class T, class Policies> bool cerge(const T& x, const interval<T, Policies>& y);
template <class T, class Policies> bool cereq(const interval<T, Policies>& x, const interval<T, Policies>& y);
template <class T, class Policies> bool cereq(const interval<T, Policies>& x, const T& y);
template <class T, class Policies> bool cereq(const T& x, const interval<T, Policies>& y);
template <class T, class Policies> bool cerne(const interval<T, Policies>& x, const interval<T, Policies>& y);
template <class T, class Policies> bool cerne(const interval<T, Policies>& x, const T& y);
template <class T, class Policies> bool cerne(const T& x, const interval<T, Policies>& y);
template <class T, class Policies> bool poslt(const interval<T, Policies>& x, const interval<T, Policies>& y);
template <class T, class Policies> bool poslt(const interval<T, Policies>& x, const T& y);
template <class T, class Policies> bool poslt(const T& x, const interval<T, Policies>& y);
template <class T, class Policies> bool posle(const interval<T, Policies>& x, const interval<T, Policies>& y);
template <class T, class Policies> bool posle(const interval<T, Policies>& x, const T& y);
template <class T, class Policies> bool posle(const T& x, const interval<T, Policies>& y);
template <class T, class Policies> bool posgt(const interval<T, Policies>& x, const interval<T, Policies>& y);
template <class T, class Policies> bool posgt(const interval<T, Policies>& x, const T& y);
template <class T, class Policies> bool posgt(const T& x, const interval<T, Policies> & y);
template <class T, class Policies> bool posge(const interval<T, Policies>& x, const interval<T, Policies>& y);
template <class T, class Policies> bool posge(const interval<T, Policies>& x, const T& y);
template <class T, class Policies> bool posge(const T& x, const interval<T, Policies>& y);
template <class T, class Policies> bool poseq(const interval<T, Policies>& x, const interval<T, Policies>& y);
template <class T, class Policies> bool poseq(const interval<T, Policies>& x, const T& y);
template <class T, class Policies> bool poseq(const T& x, const interval<T, Policies>& y);
template <class T, class Policies> bool posne(const interval<T, Policies>& x, const interval<T, Policies>& y);
template <class T, class Policies> bool posne(const interval<T, Policies>& x, const T& y);
template <class T, class Policies> bool posne(const T& x, const interval<T, Policies>& y);
/* comparison namespaces */
namespace compare {
namespace certain;
namespace possible;
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.
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. However the default comparison operators do respect this property since they throw an exception whenever [a,b] and [c,d] overlap.
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.
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.
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), rational numbers, and so
on.
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 two 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 rounding 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 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 interval library does not include I/O
operators. Printing an interval value allows a lot of customization: some
people may want to output the bounds, others may want to display the median
and the width of intervals, and so on. The example file io.cpp
shows some possibilities and may serve as a foundation in order for the user
to define her own operators.
Mixed operations with integers. When using and reusing
template codes, it is common there are operations like 2*x.
However, the library does not provide them by default because the conversion
from int to the base number type is not always correct (think
about the conversion from a 32bit integer to a single precision
floating-point number). So the functions have been put in a separate header
and the user needs to include them explicitely if she wants to benefit from
these mixed operators. Another point, there is no mixed comparison operators
due to the technical way they are defined.
Interval-aware functions. All the functions defined by
the library are obviously aware they manipulate intervals and they do it
accordingly to general interval arithmetic principles. Consequently they may
have a different behavior than the one commonly encountered with functions
not interval-aware. For example, max is defined by canonical set
extension and the result is not always one of the two arguments (if the
intervals do not overlap, then the result is one of the two intervals).
This behavior is different from std::max which returns a
reference on one of its arguments. So if the user expects a reference to be
returned, she should use std::max since it is exactly what this
function does. Please note that std::max will throw an exception
when the intervals overlap. This behavior does not predate the one described
by the C++ standard since the arguments are not "equivalent" and it allows to
have an equivalence between a <= b and &b ==
&std::max(a,b)(some particular cases may be
implementation-defined). However it is different from the one described by
SGI since it does not return the first argument even if "neither is greater
than the other".
This library was mostly inspired by previous work from Jens Maurer. Some discussions about his work are reproduced 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 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.
This material is partly based upon work supported by the National Science Foundation under NSF CAREER Grant CCR-0133599. Any opinions, findings and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the National Science Foundation (NSF).
Revised: 2004-03-11
Copyright (c) Guillaume Melquiond, Sylvain Pion, Hervé Brönnimann, 2002.
Polytechnic University.
Copyright (c) Guillaume Melquiond, 2003-2004. ENS Lyon.