Interval Computations

This chapter present the basic structures used to represent sets and the operations directly related to them.

Intervals, vectors and matrices

Ibex allows you to create interval, interval vectors, interval matrices and array of interval matrices. All these objects represent sets. These objects are respectively of the following classes:

• Interval

• IntervalVector

• IntervalMatrix

• IntervalMatrixArray.

It should be emphasized that, mainly for efficiency reasons, the previous list is not an inheritance hierarchy. So, an Interval is not a particular IntervalVector and so on, altough one might legitimately expect this (this would also be a perfectly valid design). This is in contrast with Matlab, where, basically, everything is an array of interval matrices.

Some functions return points, not intervals; this is typically the case for the function mid() that, given an interval, returns the midpoint. The returned type of Interval::mid() is theferore double. In the case of an interval vector, the function mid() returns a vector of points. In order to keep our system homogeneous, a type has been introduced to represent a vector of reals, namely, Vector and similarly for matrices or matrices arrays.

So the following types allow to represent real (floating-point) values:

• double (primitive type)

• Vector

• Matrix

• MatrixArray.

The empty set

One unpleasant consequence of our inheritance-free design is that the empty set does not have a unique representation in Ibex. The empty set is “typed” so that an empty interval is not the same object as an empty interval vector or an empty interval matrix. But, of course, they can all be interpreted as the same mathematical object, the empty set.

Worse, an empty interval vector has a size just as any other vectors (and similarly, an empty interval matrix has a certain number of rows/columns, etc.). So the empty vector of 3 components is not the same object as the empty vector of 4 components. Imposing a size for the empty vector may seem clumsy at first glance, but it actually simplifies programming in several situations. For instance, if one asks for the complementary of the empty set represented by a n-sized vector, the result will be the n-dimensional box \((-\infty,\infty)\times\ldots\times(-\infty,\infty)\). Anyway, this multiple representation is harmless in practice because we always manipulate sets of a priori known dimensions.

It holds the same with matrices and arrays of matrices.

Intervals

Intervals are represented by the class Interval (this class wraps the actual class of the underlying sub-library that implements interval arithmetic).

Here are examples of interval defintions. In particular, it is possible to give a single value to define a degenerated interval (reduced to a point). Note the usage of the constants NEG_INFINITY and POS_INFINITY for infinite bounds.

Note: In Ibex, an interval should never contain a “NaN” (not a number), whatever operation you perform with it. If this occurs, something is wrong.

	Interval x1;                      // (-oo,+oo)
	Interval x2(2);                   // [2,2]
	Interval x3(x2);                  // [2,2]
	Interval x4=x2;                   // [2,2]

	Interval x5(1,2);                 // [1,2]

	Interval x6(1,POS_INFINITY);      // [1,+oo)
	Interval x7(NEG_INFINITY,-1);     // (-oo,-1]
	Interval x8=Interval::empty_set();  // empty set

Interval Constants

Some commonly used intervals are already defined as static variables:

Interval::pi()	a thin enclosure of \(\pi\)
Interval::two_pi()	a thin enclosure of \(\pi\)
Interval::half_pi()	a thin enclosure of \(\pi/2\)
Interval::empty_set()	\(\emptyset\)
Interval::all_reals()	\((-\infty,+\infty)\)
Interval::zero()	\([0,0]\)
Interval::one()	\([1,1]\)
Interval::pos_reals()	\([0,+\infty)\)
Interval::neg_reals()	\((-\infty,0]\)

Let us print all these constants:

	output << " EMPTY_SET =\t " <<  Interval::empty_set() << endl;
	output << " PI =\t\t "      <<  Interval::pi() <<  endl;
	output << " 2 PI =\t\t "    <<  Interval::two_pi() << endl;
	output << " 1/2 PI =\t "    <<  Interval::half_pi() << endl;
	output << " ONE =\t\t "     <<  Interval::one() << endl;
	output << " ZERO =\t\t "    <<  Interval::zero() << endl;
	output << " ALL_REALS =\t " <<  Interval::all_reals() << endl;
	output << " POS_REALS =\t " <<  Interval::pos_reals() << endl;
	output << " NEG_REALS =\t " <<  Interval::neg_reals() << endl;

The output is:

 EMPTY_SET =	 [empty]
 PI =		 [3.141592653589793, 3.141592653589794]
 2 PI =		 [6.283185307179586, 6.283185307179588]
 1/2 PI =	 [1.570796326794896, 1.570796326794897]
 ONE =		 <1, 1>
 ZERO =		 <0, 0>
 ALL_REALS =	 [-inf, inf]
 POS_REALS =	 [0, inf]
 NEG_REALS =	 [-inf, 0]

Vectors or Boxes

Vectors of reals are represented by the class Vector and interval vectors (or “boxes”) are represented by the class IntervalVector.

Creation and initialization of interval vectors can either be done separately (you first create a vector of the desired size and then fill it with the intervals) or simultaneously, by giving an array of double in argument of the constructor. The following example creates twice ([0,1],[2,3],[4,5]):

	// creates ([0,1],[2,3],[4,5]) in several steps.
	IntervalVector x1(3); // so far x1 is ((-oo,oo),(-oo,oo),(-oo,oo))
	x1[0]=Interval(0,1);  // we initialize each component
	x1[1]=Interval(2,3);
	x1[2]=Interval(4,5);

	// creates the same vector in two steps only
	double _x2[3][2]={{0,1},{2,3},{4,5}};
	IntervalVector x2(3,_x2);

This is the same for vectors:

	// creates (0.1,0.2,0.3) in several steps.
	Vector x1(3);  // so far x1 is (0,0,0)
	x1[0]=0.1;     // we initialize each component
	x1[1]=0.2;
	x1[2]=0.3;

	// creates the same vector in two steps only
	double _x2[3]={0.1,0.2,0.3};
	Vector x2(3,_x2);

You can also create an (interval) vector by duplicating n times the same double/interval and copy another (interval) vector:

	IntervalVector x2(3,Interval(1,2));         // create ([1,2],[1,2],[1,2])
	IntervalVector x3(3,0.1);                   // create ([0.1,0.1],[0.1,0.1],[0.1,0.1])
	Vector         x4(3,0.1);                   // create (0.1,0.1,0.1)

	IntervalVector x5(x2);                      // create a copy of x2
	Vector x6(x4);                              // create a copy of x4

Finally, you can create an empty interval vector of a given size (see The empty set) or create a degenerated box from a vector of reals:

	IntervalVector x1=IntervalVector::empty(3); // create a vector of 3 empty intervals

	Vector x(3,0.1);
	IntervalVector x2(x);                       // create ([0.1,0.1],[0.1,0.1],[0.1,0.1])

Matrices and Array of matrices

They work exactly the same way as vectors. Here are first listed different ways to build a Matrix:

	// create [(1,2,3);(4,5,6)] in several steps.
	Matrix m1(2,3);   // a 2x3 matrix filled with zeros
	m1[0][0]=1;       // we initialize each component
	m1[0][1]=2;
	m1[0][2]=3;
	m1[1][0]=4;
	m1[1][1]=5;
	m1[1][2]=6;
	output << "m1=" << m1 << endl << endl;

	// create the same matrix in two steps only
	double _m2[2*3]={1,2,3,
			 4,5,6};
	Matrix m2(2,3,_m2);
	output << "m2=" << m2 << endl << endl;

	// create a 2x3 matrix filled with ones
	Matrix m3(2,3,1.0);
	output << "m3=" << m3 << endl << endl;

	// the same matrix using a built-in static function
	Matrix m4=Matrix::ones(2,3);

	// create a copy of m3
	Matrix m5(m3);
	output << "m5=" << m5 << endl;

The output is:

m1=((1 ; 2 ; 3)
(4 ; 5 ; 6))

m2=((1 ; 2 ; 3)
(4 ; 5 ; 6))

m3=((1 ; 1 ; 1)
(1 ; 1 ; 1))

m5=((1 ; 1 ; 1)
(1 ; 1 ; 1))

Now, for interval matrices, we shall take as example a matrix that frequently appears, the “perturbed” idendity, that is:

\[\begin{split}\begin{pmatrix} [1-\varepsilon,1+\varepsilon] & \ldots & [-\varepsilon,\varepsilon] \\ \vdots & \ddots & \vdots\\ [-\varepsilon,\varepsilon] & \ldots & [1-\varepsilon,1+\varepsilon] \end{pmatrix}\end{split}\]

	double eps=1e-02;

	// create in several steps.
	IntervalMatrix m1(3,3);                   // a 3x3 matrix filled with (-oo,oo)
	m1=Matrix::eye(3);                        // set m1 to the identity matrix
	m1+=Interval(-eps,eps)*Matrix::ones(3);   // add [-eps,eps] to each component
	output << "m1=" << m1 << endl << endl;

	// create the same matrix with a matrix of double
	// warning: the matrix of double has a number of rows
	// equal to the total number of components (9) and
	// 2 columns (left bound, right bound).
	double _m2[3*3][2]={{1-eps,1+eps},{-eps,eps},{-eps,eps},
	                   {-eps,eps},{1-eps,1+eps},{-eps,eps},
	                   {-eps,eps},{-eps,eps},{1-eps,1+eps}};

	IntervalMatrix m2(3,3,_m2);
	output << "m2=" << m2 << endl << endl;

	// create a 3x3 matrix filled with [-eps,eps]
	IntervalMatrix m3(3,3,Interval(-eps,eps));
	m3+=Matrix::eye(3);
	output << "m3=" << m3 << endl << endl;

	// create a copy of m3
	IntervalMatrix m4(m3);
	output << "m4=" << m4 << endl;

The output is:

m1=(([0.9899999999999999, 1.010000000000001] ; [-0.01000000000000001, 0.01000000000000001] ; [-0.01000000000000001, 0.01000000000000001])
([-0.01000000000000001, 0.01000000000000001] ; [0.9899999999999999, 1.010000000000001] ; [-0.01000000000000001, 0.01000000000000001])
([-0.01000000000000001, 0.01000000000000001] ; [-0.01000000000000001, 0.01000000000000001] ; [0.9899999999999999, 1.010000000000001]))

m2=(([0.9899999999999999, 1.010000000000001] ; [-0.01000000000000001, 0.01000000000000001] ; [-0.01000000000000001, 0.01000000000000001])
([-0.01000000000000001, 0.01000000000000001] ; [0.9899999999999999, 1.010000000000001] ; [-0.01000000000000001, 0.01000000000000001])
([-0.01000000000000001, 0.01000000000000001] ; [-0.01000000000000001, 0.01000000000000001] ; [0.9899999999999999, 1.010000000000001]))

m3=(([0.9899999999999999, 1.010000000000001] ; [-0.01000000000000001, 0.01000000000000001] ; [-0.01000000000000001, 0.01000000000000001])
([-0.01000000000000001, 0.01000000000000001] ; [0.9899999999999999, 1.010000000000001] ; [-0.01000000000000001, 0.01000000000000001])
([-0.01000000000000001, 0.01000000000000001] ; [-0.01000000000000001, 0.01000000000000001] ; [0.9899999999999999, 1.010000000000001]))

m4=(([0.9899999999999999, 1.010000000000001] ; [-0.01000000000000001, 0.01000000000000001] ; [-0.01000000000000001, 0.01000000000000001])
([-0.01000000000000001, 0.01000000000000001] ; [0.9899999999999999, 1.010000000000001] ; [-0.01000000000000001, 0.01000000000000001])
([-0.01000000000000001, 0.01000000000000001] ; [-0.01000000000000001, 0.01000000000000001] ; [0.9899999999999999, 1.010000000000001]))

Container operations

This section presents operations on interval vectors and interval matrices viewed as containers.

Most of the operations are also possible with vectors (Vector) and matrices (Matrix) by replacing Interval with double and so on. They are not detailed again here, please refer to the API.

Container operations for an interval vector x

Return type	C++ code	Meaning
int	x.size()	The size (number of components)
Interval&	x[i] (int i)	The ith component. By reference.
IntervalVector	x.subvector(i,j) (int i, int j)	Return the subvector of length (j-i+1) starting at index i and ending at index j (both included). By copy. [x] must be non-empty.
void	x.resize(n) (int n)	Resize the vector. If the vector has been enlarged the extra components are set to \((-\infty,+\infty)\) even if the vector is empty.
void	x.put(i, y) (int i, IntervalVector& y)	Write the vector [y] in [x] at index i.
void	x.clear()	Set all the elements to [0,0] (even if [x] is empty). Emptiness is “overridden”.
void	x.init(y) (Interval& y)	Set all the elements to [y] (even if [x] is empty). Emptiness is “overridden”.

Container operations for an interval matrix m

Return type	C++ code	Meaning
int	m.nb_rows()	Number of rows
int	m.nb_cols()	Number of column
Interval&	m[i][j] (int i, int j)	The (i,j)th entry. By reference.
IntervalVector&	m.row(i) (int i)	The ith row. By reference (fast).
IntervalVector	m.col(j) (int j)	The jth column. By copy (slow, in O(n))
IntervalMatrix	m.rows(i1,i2) (int i1, int i2)	The submatrix obtained by selecting rows between index i1 and i2. By copy.
IntervalMatrix	m.cols(j1,j2) (int j1, int j2)	The submatrix obtained by selecting columns between index j1 and j2. By copy.
IntervalMatrix	m.submatrix(i1,i2,j1,j2) (int i1, int i2, int j1, int j2)	The submatrix by selecting rows from index i1 to i2 (both included) and columns from index j1 to j2 (both included).
void	m.resize(n1,n2) (int n1, int n2)	Resize the matrix to a n1xn2 matrix. Extra components are set to \((-\infty,+\infty)\) even if the matrix is empty.
void	m.set_row(i, y) (int i, IntervalVector& y)	Write the vector [y] in the ith row of [m].
void	m.set_col(i, y) (int i, IntervalVector& y)	Write the vector [y] in the ith column of [m].
void	m.put(i,j,m2) (int i, int j, IntervalMatrix& m2)	Write the matrix [m2] in [m] at (i,j) and subsequent indices
void	m.put(i,j,v,b) (int i, int j, IntervalVector& v, bool b)	Write the vector at (i,j) either in row (if b==true) or in column (if b==false).
void	m.clear()	Set all the elements to [0,0] (even if [m] is empty). Emptiness is “overridden”.
void	m.init(y) (Interval& y)	Set all the elements to [y] (even if [m] is empty). Emptiness is “overridden”.

Set-membership operations

The operations described here are valid for any “sets” [x] and [y] and any “point” p, where a “set” designates here one object of either of the following classes:

• Interval

• IntervalVector

• IntervalMatrix

• IntervalMatrixArray.

Similarly, a “point” either refers to a double or an instance of:

• Vector

• Matrix

• MatrixArray.

Given equidimensional sets [x], [y] and a point p, the following table gives the set-membership predicates currently available in Ibex (from release 2.8) and their mathematical intepretation:

C++ code	Meaning	available from
x==y	\([x]=[y]\)
x!=y	\([x]\neq [y]\)
x.is_empty()	\([x]=\emptyset\)
x.is_subset(y)	\([x]\subseteq [y]\)
x.is_strict_subset(y)	\([x]\subseteq [y]\wedge [x]\neq [y]\)
x.is_interior_subset(y)	\([x]\subseteq \mathring{[y]}\)	release 2.8
x.is_strict_interior_subset(y)	\([x]\subseteq \mathring{[y]}\wedge [x]\neq [y]\)	release 2.8
x.is_superset(y)	\([x]\supseteq [y]\)
x.is_strict_superset(y)	\([x]\supseteq [y]\wedge [x]\neq [y]\)
x.contains(p)	\(p\in [x]\)
x.interior_contains(p)	\(p\in \mathring{[x]}\)	release 2.8
x.intersects(y)	\([x]\cap [y]\neq\emptyset\)	release 2.8
x.overlaps(y)	\(\mathring{[x]}\cap \mathring{[y]}\neq\emptyset\)	release 2.8
x.is_disjoint(y)	\([x]\cap [y]=\emptyset\)	release 2.8

Possible set-membership operations between sets are:

C++ code	Meaning
x&y	\([x]\cap [y]\)
x|y	\(\square([x]\cup [y])\)
x.set_empty()	\([x]\leftarrow \emptyset\)
x=y	\([x]\leftarrow [y]\)
x&=y	\([x]\leftarrow ([x]\cap [y])\)
x|=y	\([x]\leftarrow \square([x] \cup [y])\)

Finally, given two interval vectors [x] and [y] one can call:

cart_prod(x,y)

to create the Cartesian product.

Complementary and set difference

These set operations are only available for intervals (from release 2.8) and interval vectors.

In the case of intervals, the function complementary stores the complementary of *this in two intervals c1 and c2. More exactly, the result is the closure of the complementary since open bounds are not representable (except for infinity and minus infinity). E.g, the “complementary” of \((-\infty,0]\) is \([0,+\infty)\).

	Interval x(0,1);
	Interval c1,c2; // to store the result

	int n=x.complementary(c1,c2);
	output << "complementary of " << x << " = " << c1;
	if (n>1) output << " and " << c2;
	output << endl;

The output is:

complementary of [0, 1] = [-inf, 0] and [1, inf]

In the case of interval vectors, the function complementary calculates the complementary under the form of a union of non-overlapping interval vectors, and stores this union into an array that is allocated by the function itself (the variable result in the code below). The function returns the size of the union/array.

To illustrate this, let us first build a function that prints the complementary of an interval vector:

void print_compl(const IntervalVector& x) {
	IntervalVector* result;
	int n=x.complementary(result);
	output << "complementary of " << x << " = " << endl;
	for (int i=0; i<n; i++) {
	output << "\t" << result[i] << endl;
	}

	delete[] result; // don't forget to free memory!
}

We can call it now with different vectors. Note that when the vector is the empty set with n components, the complementary is a n-dimensional box: \((-\infty,\infty)\times\ldots\times(-\infty,\infty)\). Note also that if the difference is empty, result is an array of one element set to the empty box. It is not a zero-sized array containing no element (this is illegal in ISO C++). However the returned number is 0 (not 1). The interesting point is that you can call delete[] safely, in all cases.

	print_compl(IntervalVector::empty(3));

	print_compl(IntervalVector(3));

	print_compl(IntervalVector(3,Interval(0,1)));

The output is:

complementary of empty vector = 
	([-inf, inf] ; [-inf, inf] ; [-inf, inf])
complementary of ([-inf, inf] ; [-inf, inf] ; [-inf, inf]) = 
complementary of ([0, 1] ; [0, 1] ; [0, 1]) = 
	([-inf, 0] ; [-inf, inf] ; [-inf, inf])
	([1, inf] ; [-inf, inf] ; [-inf, inf])
	([0, 1] ; [-inf, 0] ; [-inf, inf])
	([0, 1] ; [1, inf] ; [-inf, inf])
	([0, 1] ; [0, 1] ; [-inf, 0])
	([0, 1] ; [0, 1] ; [1, inf])

The set difference works exactly the same way except that the function takes as first argument another set [y].

Example with intervals:

	// set difference between two intervals
	Interval x(0,3);
	Interval y(1,2);
	Interval c1,c2; // to store the result

	int n=x.diff(y,c1,c2);
	output << x << " \\ " << y << " = " << c1;
	if (n>1) output << " and " << c2;
	output << endl;

[0, 3] \ [1, 2] = [0, 1] and [2, 3]

Example with interval vectors:

	// set difference between two boxes
	IntervalVector x(2,Interval(0,3));
	IntervalVector y(2,Interval(1,2));
	IntervalVector* result; // to store the result

	int n=x.diff(y,result);
	output << x << " \\ " << y << " = " << endl;
	for (int i=0; i<n; i++) {
		output << "\t" << result[i] << endl;
	}
	delete[] result; // don't forget to free memory!

([0, 3] ; [0, 3]) \ ([1, 2] ; [1, 2]) = 
	([0, 1] ; [0, 3])
	([2, 3] ; [0, 3])
	([1, 2] ; [0, 1])
	([1, 2] ; [2, 3])

Geometric operations

Here are the functions that can be applied to an interval [x], seen as a segment of the line.

Return type	C++ code	Meaning
double	x.lb()	\(\underline{x}\), the lower (left) bound of [x]
double	x.ub()	\(\overline{x}\), the upper (right) bound of [x]
double	x.diam()	diameter, \(|\overline{x}-\underline{x}|\)
double	x.rad()	radius, half of the diameter
double	x.mid()	the midpoint, (\((\underline{x}+\overline{x})/2\))
Interval	x.inflate(eps) (double eps)	an interval with the same midpoint and radius increased by eps
double	distance(x,y) (Interval& y)	the (Hausdorff) distance between [x] and [y]. The distance is 0 iff [x]==[y]. Otherwise, it returns the minimal value by which one of the intervals ([x] or [y)) has to be inflated so that it entirely overlaps the other interval.
double	x.rel_distance(y) (Interval& y)	the “relative” distance, that is, distance(x,y)/x.diam()
bool	x.is_unbounded()	true iff [x] has one of its bounds infinite.
bool	x.is_bisectable()	true iff [x], as an interval of floating point numbers, can be bisected in two non-degenerated intervals of floating point numbers. The empty interval or an interval of two consecutive floating points are not bisectable.

Here are the functions that can be applied to an interval vector [x], seen as a box.

Return type	C++ code	Meaning
Vector	x.lb()	lower-left corner (vector of lower bounds of [x])
Vector	x.ub()	upper-right corner (vector of upper bounds of [x])
Vector	x.diam()	vector of diameters, \(|\overline{x_i}-\underline{x_i}|\)
double	x.min_diam()	minimal diameter, among all components of [x]
double	x.max_diam()	maximal diameter, among all components of [x]
int	x.extr_diam_index(b) (bool b)	the index of a component with minimal (if b==true) or maximal (if b==false) diameter.
void	x.sort_indices(b, tab) (bool b, int tab[])	Write into tab the indices of all the components, sorted by increasing (if b==true) /decreasing (if b==false) diameters. The array tab must have been allocated before calling this function.
Vector	x.rad()	vector of radii (halves of diameters)
Vector	x.mid()	the midpoint, (\((\underline{x}+\overline{x})/2\))
double	x.volume()	the volume of the box
double	x.perimeter()	the perimeter of the box
bool	x.is_flat()	true if the volume is null (one dimension is degenerated)
IntervalVector	x.inflate(eps) (double eps)	a box with the same midpoint and each radius increased by eps%
double	distance(x,y) (IntervalVector& y)	the (Hausdorff) distance between [x] and [y]. The distance is 0 iff [x]==[y]. Otherwise, it returns the minimal value by which one of the boxes ([x] or [y]) has to be inflated so that it entirely overlaps the other box.
double	x.rel_distance(y) (IntervalVector& y)	the “relative” distance, that is, distance(x,y)/x.max_diam()
bool	x.is_unbounded()	true iff [x] has one of its bounds infinite.
bool	x.is_bisectable()	true iff at least one component of [x] is bisectable (see above).

Bisection

Bisecting a box is a fundamental operation in Ibex. The choice of the component to be bisected is often critical and defining a strategy for chosing this component is the purpose of Bisectors.

To bisect the ith component, the bisect function can be used. The first arugment is the index of the dimension i along which the box has to be splitted. The second argument is optionnal and gives at which point the ith interval is bisected. This argument is entered as a ratio of the interval diameter, 0.5 means “midpoint” (and the default value is precisely 0.5). Example:

	IntervalVector x(3,Interval(0,1)); // [0,1]x[0,1]x[0,1]

	std::pair<IntervalVector,IntervalVector> p = x.bisect(1,0.4); // bisect the second component with ratio 0.4

	output << "first box=" << p.first << endl;
	output << "second box=" << p.second << endl;

The output is:

first box=([0, 1] ; [0, 0.4000000000000001] ; [0, 1])
second box=([0, 1] ; [0.4, 1] ; [0, 1])

Miscellaneous

Given an interval [x]:

C++ code	Meaning
x.mig()	Mignitude, \(\min_{x\in[x]} |x|\) (a double). Also exists for interval vectors and matrices (gives a vector or matrix of mignitudes).
x.mag()	Magnitude, \(\max_{x\in[x]} |x|\) (a double). Also exists for interval vectors and matrices (gives a vector or matrix of magnitudes).
integer(x)	The largest interval [a,b] of integers included in [x].

Finally, for an interval vector [v]:

v.random()

A random point inside x

Interval arithmetic

Interval arithmetic is the main device upon wich Ibex is built.

For a complete introduction of interval arithmetic and interval analysis, see [Moore 1966], [Neumaier 1990] or [Jaulin et al. 2001].

The interval arithmetic defines for each elementary function \(f:\mathbb{R}^n\to\mathbb{R}\) an interval function \([f]:\mathbb{IR}^n\to\mathbb{IR}\) according to the formula:

\[[f]([x]) \supseteq \square\{f(\alpha), \ \alpha\in [x] \}\]

where \(\square\) (the “hull” symbol) in front of a set means “the smallest box enclosing this set”. The hull symbol is here to highlight the fact that an enclosure of the possibly discontinuous set is returned. Note that, among all the functions below, only tan and sign fail to be continuous. So, in the other cases, the “hull” operator can be ommitted. The possible roundoff error accounts for the \(\supseteq\) symbol.

The function [f] is called an “interval extension” of the original function f (applied on real numbers).

We start by detailing the binary operators (+,-,*,/) and then proceed with the nonlinear ones.

Addition, subtraction, multiplication

The following tables sumarize the basic operations of linear algebra in Ibex. It is important to notice that (so far?), theses operations are not optimized, as compared to specialized softwares like Intlab. They are all based on naive implementations.

• Symbols x and y denote scalars (double or instances of Interval);

• Symbols v and w denote vectors (instances of Vector or IntervalVector)

• Symbols A and B denotes matrices (instances of Matrix or IntervalMatrix).

Notice that it is not possible (at least in the current release) to appply such operations with objects of type IntervalMatrixArray, the latter being considered as a pure “set” or, more pragmmatically, as a matrix container. In other words, arrays of interval matrices are not a 3rd order tensor.

Possible additions and subtractions (with optionally assignments) are:

-x	x+y	x-y	x+=y	x-=y
-v	v+w	v-w	v+=w	v-=w
-A	A+B	A-B	A+=B	A-=B

Possible multiplications are:

C++ code	Meaning
x*y	multiplication of two “scalars” (double/intervals)
x*=y	multiplication and assignment
x*v	scalar multiplication of a vector
x*A	scalar multiplication of a matrix
v*w	dot product (\(v^Tw\))
hadamard_product(v,w)	\((v_1w_1,\ldots,v_nw_n)\)
outer_product(v,w)	outer product (\(vw^T\))
A*x	matrix-vector product
A*B	matrix product

Division

First, it is possile to use division between two interval/doubles using C++ operator overloading, as for the previous operators. You can either write x/y (division) or x/=y (division and assignment).

However, as a set-membership operator, and because the division of reals is non-continuous, dividing two intevals may result in a union of two intervals. E.g.,

\[\,\![2,3] / [-1,2] = \{ x/y, x\in[2,3], y\in[-1,2] \} = (-\infty,-2]\cup[1,+\infty)\]

This variant of the interval division is called “generalized division” in the litterature. The previous division operator / simply gives (in Ibex) the hull of this union. So, in C++, [2,3] / [-1,2] gives \((-\infty,\infty)\). If you want to handle the possible outcome of two intervals, you have to call the “div2” function.

Given four intervals [x],[y], [out1] and [out2]:

div2(x,y,out1,out2)

will store the union representing [x]/[y] in [out1] and [out2]. If the result is a single interval, [out2] is set to the empty interval. If the result is the empty set, both [out1] and [out2] are set to the empty interval.

Note: Contrary to the “cset” theory, the result is empty if y=[0,0] (whatever [x] is).

For convenience, there is also a function in the Interval class that simultaneously performs divsion and intersection:

bool div2_inter(const Interval& x, const Interval& y, Interval& out2)

This function sets *this to the intersection of itself with the division of two others, [x] and [y].

In return, *this and [out2] contains the lower and upper part respectively of the division. If the result of the generalized division and intersection is a single interval, [out2] is set to the empty interval.

The function returns true if the intersection is non empty.

Example:

	Interval intv(-10,10);
	Interval out2;
	bool result=intv.div2_inter(Interval(2,3), Interval(-1,2), out2);
	output << "the intersection is " << (result? "not":"") << " empty" << endl;
	output << "left part=" << intv << " right part=" << out2 << endl;

the intersection is not empty
left part=[-10, -2] right part=[1, 10]

Non-linear elementary functions

The following operations are allowed for an interval [x].

Power and roots

C++ code	Meaning
sqr(x)	\([x]^2\)
sqrt(x)	\(\sqrt{[x]}\)
pow(x,n)	\([x]^n\)
pow(x,y)	\([x]^{[y]} = e^{[y]\log([x])}\)
root(x,n)	\(\sqrt[n]{[x]}\)

Exponenial, logarithm

exp(x)

log(x)

Trigonometric and hyperbolic functions

cos(x)	sin(x)	tan(x)
acos(x)	asin(x)	atan(x)
cosh(x)	sinh(x)	tanh(x)
acosh(x)	asinh(x)	atanh(x)
atan2(y,x)

Max, min and miscellaneous

The sign function below is the interval extension of

\[\begin{split}sign(\alpha) := \left\{\begin{array}{ll} -1 & \mbox{if} \ \alpha<0 \\ 0 & \mbox{if} \ \alpha=0 \\ 1 & \mbox{if} \ \alpha>0 \\ \end{array}\right.\end{split}\]

C++ code	Meaning
abs(x)	interval extension of the absolute value. Not to be confused with magnitude; e.g. with [x]=[-2,1], abs([x])=[0,2] and the magnitude is 2.
max(x,y)	interval extension of the maximum. Not to be confused with \(\max ([x]\cup [y])\); e.g., max([0,3],[1,2])=[1,3].
min(x,y)	interval extension of the minimum. Not to be confused with \(\min ([x]\cup [y])\); e.g., min([0,3],[1,2])=[0,2].
sign(x)	(see above)

Backward arithmetic

The “backward” or “relationnal” arithmetic consists in contraction operators for elementary constraints like

\[y= x_1+x_2.\]

More precisely, given a function \(f:\mathbb{R}\to\mathbb{R}\) (\(\sin, \cos, \exp\, \ldots\)) the backward operator of f is an interval function \([f]_{bwd}:\mathbb{IR}^2\to\mathbb{IR}\) that satisfies:

\[[f]_{bwd}([x],[y]) \supseteq \{ x\in [x] \ | \ \exists y\in [y], \quad y=f(x) \}.\]

Similarly, given a binary function \(f:\mathbb{R}^2\to\mathbb{R}\) (\(+, -, \times, \max, \ldots\)) the backward operator is a function \([f]_{bwd}:\mathbb{IR}^2\to\mathbb{IR}\) that satisfies:

\[[f]_{bwd}([x]_1,[x]_2,[y])\supseteq\{ (x_1,x_2)\in [x]_1\times[x]_2 \ | \ \exists y\in[y],\ y=f(x_1,x_2)\}.\]

In practice, backward operators in Ibex do not return an interval, they directly contracts the input interval [x] (in case of unary functions), that is:

\[[x] \leftarrow [f]_{bwd}([x],[y])\]

or the two intervals [x1], [x2] (in case of binary functions):

\[([x]_1,[x]_2) \leftarrow [f]_{bwd}([x]_1,[x]_2,[y]).\]

So the bwd_* functions below all return void.

Note: One important feature in Ibex is the ability to contract with respect to any constraints and, in particular, with constraints under the form

\[y=f(x_1,\ldots,x_n)\]

where f is an arbirary function. So, as a user, there is probably little interset for you to call the low-level routines presented here. Moreover, these routines only contracts the “pre-image argument” [x]; they don’t contract the “image argument” [y]. This is in contrast with high-level contractors like Forward-Backward that do both contractions (on [x] and [y]).

C++ code	Relation
bwd_add(y,x1,x2)	\(y=x_1+x_2\)
bwd_sub(y,x1,x2)	\(y=x_1-x_2\)
bwd_mul(y,x1,x2)	\(y=x_1\times x_2\)
bwd_div(y,x1,x2)	\(y=x_1/x_2\)
bwd_sqr(y,x)	\(y=x^2\)
bwd_sqrt(y,x)	\(y=\sqrt{x}\)
bwd_pow(y,n,x)	\(y=x^n\)
bwd_pow(y,x1,x2)	\(y=x_1^{x_2}\)
bwd_root(y,n,x)	\(y=\sqrt[n]{x}\)
bwd_exp(y,x)	\(y=\exp{x}\)
bwd_log(y,x)	\(y=\ln{x}\)
bwd_cos(y,x)	\(y=\cos{x}\)
bwd_sin(y,x)	\(y=\sin{x}\)
bwd_tan(y,x)	\(y=\tan{x}\)
bwd_acos(y,x)	\(y=\arccos{x}\)
bwd_asin(y,x)	\(y=\arcsin{x}\)
bwd_atan(y,x)	\(y=\arctan{x}\)
bwd_cosh(y,x)	\(y=\cosh{x}\)
bwd_sinh(y,x)	\(y=\sinh{x}\)
bwd_tanh(y,x)	\(y=\tanh{x}\)
bwd_acosh(y,x)	\(y=\operatorname{arcosh}(x)\)
bwd_asinh(y,x)	\(y=\operatorname{arsinh}(x)\)
bwd_atanh(y,x)	\(y=\operatorname{artanh}(x)\)
bwd_atan2(y,x1,x2)	\(y=atan2(x_1,x_2)\)
bwd_abs(y,x)	\(y=|x|\)
bwd_sign(y,x)	\(y=sign(x)\) (see def. above)
bwd_integer(y,x)	\(y=x \wedge x\in\mathbb{N}\)
bwd_min(y,x1,x2)	\(y=\min(x_1,x_2)\)
bwd_max(y,x1,x2)	\(y=\max(x_1,x_2)\)

Inner arithmetic

The inner arithmetic [Chabert & Beldiceanu 2010] [Araya et al 2014] consists in two types of operators. Given a function \(f:\mathbb{R}^n\to\mathbb{R}\):

• the forward inner operator is an interval function \(]f[:\mathbb{IR}^n\to\mathbb{IR}\) that satisfies the formula:

  \[]f[([x]) \subseteq \{f(\alpha), \ \alpha\in [x] \}\]

• the backward inner operator is an interval function \(]f[_{bwd}:\mathbb{IR}^{n+1}\to\mathbb{IR}\) that satisfies the formula:

  \[]f[_{bwd}([x],[y]) \subseteq \{ x\in [x] \ | \ \exists y\in [y], \quad y=f(x) \}.\]

Notice that the inclusion symbol has just been reversed, as compared to the classical interval (forward) and backward arithmetic.

Forward operators

Note: The inner operators are not implemented for all elementary functions yet.

C++ code	Function
iadd(x1,x2)	x1+x2
isub(x1,x2)	x1-x2
imul(x1,x2)	x1*x2
idiv(x1,x2)	x1/x2
isqr(x)	x^2
iminus(x)	-x
ilog(x)	ln(x)
iexp(x)	exp(x)
iacos(x)	cos(x)
iasin(x)	sin(x)
iatan(x)	tan(x)

Backward operators

Let us denote

\[\mathcal{S}:=\{ x\in [x] \ | \ \exists y\in [y], \quad y=f(x) \}.\]

The backward operator also takes an optional argument xin that represents an “inner” box, that is a box which already satisfies

\[[x]_{in} \subseteq \mathcal{S}.\]

If this argument is set, the operator guarantees that the contracted box [x] will enclose [xin]. The operator works in this case as an inflator, the inner box [xin] is indeed “inflated” to a larger inner box [x]. More precisely, we have:

\[[x]_{in} \ \subseteq \ ]f[_{bwd}([x],[y],[x]_{in}) \subseteq ([x] \cap \mathcal{S}).\]

Note: The inner operators are not implemented for all elementary functions yet.

C++ code		Relation
ibwd_add(y,x1,x2,xin1,xin2)		y=x1+x2
ibwd_sub(y,x1,x2,xin1,xin2)		y=x1-x2
ibwd_mul(y,x1,x2,xin1,xin2)		y=x1*x2
ibwd_div(y,x1,x2,xin1,xin2)		y=x1/x2
ibwd_sqr(y,x,xin)		y=x^2
ibwd_sqrt(y,x,xin)		\(y=\sqrt{x}\)
ibwd_minus(y,x,xin)		y=-x
ibwd_abs(y,x,xin)		\(y=|x|\)
ibwd_pow(y,x,p,xin)		y=x^p
ibwd_log(y,x,xin)		y=ln(x)
ibwd_exp(y,x,xin)		y=exp(x)
ibwd_cos(y,x,xin)		y=cos(x)
ibwd_sin(y,x,xin)		y=sin(x)
ibwd_tan(y,x,xin)		y=tan(x)