Final update: March 29, 2016

# Interval Arithmetic Library

Masahide Kashiwagi

## 2. Features

• This library provides rigorous arithmetics for closed intervals represented by its upper and lower bounds as boost.interval.
• The upper and lower bounds are not restricted into the double type.
• boost.interval employs two members in its interval class as interval<T,P>; T represents the type of the pair of bounds, and P (called as the policy) indicates several specifications, e.g., of a rounding mode and so on. This is an "advanced" data structure, however, the simple data type as interval<T> is employed in KV library abolishing the member of policy because of several technical reasons.
• The usage of almost all functions in KV library follows boost.interval.
• KV library supports the interconversion between binary and decimal numbers with changing the current rounding mode appropriately (whereas boost.interval does not support such interconversion). Owing to this function, for example, cout << x; displays the lower and upper bound of x which are rounded down and rounded up, respectively. Moreover, when writing x="0.1";, an interval strictly including 0.1 is substituted into x.
• KV library supports the interval arithmetics for several mathematical functions, e.g., of exp, log, sin, cos, tan, sinh, cosh, tanh, asin, acos, atan, asinh, acosh, atanh, expm1, log1p, abs, and pow, whereas boost.interval does not support such mathematical functions without combining with other libraries (e.g., crlibm).
• A High-precision interval arithmetic with about 32 significant figures can be performed by assigning the member T (in interval<T>) to the dd type. All functions (including the interconversion between binary and decimal numbers and the above mathematical functions) accurately works for the dd type with a suitable precision.

## 3. Supported environment

• A computational environment where C++ program and Boost can work, is required.
• In order to obtain correct numerical results using rdouble.h (to be introduced later), one of the following conditions is required:
• "fesetround" given by fenv.h works. Recent gcc Compilers that is C99 standard, supports "fesetround".
• "_controlfp" given by float.h supported by Microsoft Visual C++ works.
Note that the option of -DKV_FASTROUND is available to perform faster interval arithmetics when using an Intel's CPU and a 64bit operating system. See hwround.h for details. A Unix operating system (64bit) is recommended for the use of KV library.

## 4. Files

interval.hpp
(requiring convert.hpp and constants.hpp)
rdouble.hpp
(requiring conv-double.hpp, hwround.hpp, rdouble-hwround.hpp, and rdouble-nohwround.hpp)

It is required to include rdouble.hpp together with interval.hpp for the use of interval arithmetic with the double precision. The computation using interval.hpp without including rdouble.hpp does not change rounding modes in a computation process, and therefore returns unverified results.

## 5. How to use

### 5.1 Basic instruction

Inclusion
```#include <kv/interval.hpp>
```
enables us to use interval arithmetics. Note that, if the type of the bounds of an interval requires the change of a rounding mode (e.g., the double type or the dd type), an additional header file is required to be included. For example, when the value type in an interval is double, one has to additionally include kv/rdouble.hpp as
```#include <kv/interval.hpp>
#include <kv/rdouble.hpp>
```

The type of interval in the namespace of kv is defined with a template style. For example, interval<double> represents the interval, the value type in which is double. Therefore, the variable declaration for the interval can be done as

```kv::interval<double> x;
```
"kv::" can be omitted by declaring "using namespace kv;". Otherwise, after doing "typedef" as
```typedef kv::interval<double> itv;
```
declaring a variable as
```itv x;
```
is convenient for a simple programing. Hereafter, this notation is used for simplicity

### 5.2 Declaration and substitution

The usual variable declaration is done as
```    itv x;
```
as mentioned above. Variable declarations with initializing data are done as
```    itv x(1.);
itv x(1., 2.);
```
An error occurs by the form
```    itv x = 1.;
```
because the constructor of one argument is explicit. No errors occur by the forms
```    itv x = (itv)1.;
itv x = itv(1.);
itv x = itv(1., 2.);
```

The assignment operator for an interval type is overloaded, therefore no errors occur by the forms

```    itv x;
x = 1.;
x = (itv)1.;
x = itv(1.);
x = itv(1., 2.);
```
Remark.   The type of input arguments in initialization and substitution for an interval (e.g., the above "1." or "2.") can be assigned to a type that is convertible to the type in the objective interval. For example, "x = 1." works for the dd-interval type (interval<dd>). See here (in Japanese) for details.
To be compatible with boost.interval, the form like
```    x.assign(1., 2.);
```
works.

The initialization and the substitution of characters are supported as

```    itv x("0.1");
itv x("0.1", "0.2");
```
or
```    itv x;
x = "0.1";
x = (itv)"0.1";
x = itv("0.1");
x = itv("0.1", "0.2");
```
In this case, when writing x = "0.1";, the lower bound of x is the maximum double-number less than or equal to 0.1, and the upper bound of x is the minimum double-number more than or equal to 0.1; therefore x strictly include 0.1. Likewise, ("0.1", "0.2") denotes the interval whose lower bound is maximum double-number less than or equal to 0.1, and whose upper bound is the minimum value more than or equal to 0.2.

Note that such translations of characters are not fast. Therefore, it is recommended to prevent any loops (be it, for, while, and so on.) from having such translations.

### 5.3 Display

```    std::cout << z << "\n";
```
works for the display of an interval z, where the binary‐to‐decimal conversion is done with an appropriate rounding for the correctness of the display. "cout.precision" is available as
```    std::cout.precision(17);
std::cout << z << "\n";
```
which display z with the 17-digits precision.

### 5.4 Four arithmetic operations

The following codes work:
```    itv x, y, z;
x = 1.; y = 2.;

z = x + y;
z = x - y;
z = - x;
z = x * y;
z = x / y;

z = x + 1.;
z = 2. + x;
z = x + "0.1";

x += y;
x += 1.;
x += "0.1";
```
The type of constants in an operation with an interval (e.g., the above "1.", "2.", and "0.1") can be assigned to a type convertible to the type in the objective interval.

### 5.5 Mathematical functions

The following codes work:
```    std::cout << sqrt(itv(2.5, 3.5)) << "\n";
std::cout << exp(itv(2.5, 3.5)) << "\n";
std::cout << expm1(itv(-0.25, 0.25)) << "\n";
std::cout << log(itv(0.75, 1.25)) << "\n";
std::cout << log(itv(0., 1.)) << "\n";
std::cout << log1p(itv(-0.25, 0.25)) << "\n";
std::cout << sin(itv(-0.25, 0.25)) << "\n";
std::cout << cos(itv(-0.25, 0.25)) << "\n";
std::cout << tan(itv(-0.25, 0.25)) << "\n";
std::cout << atan(itv(-0.25, 0.25)) << "\n";
std::cout << asin(itv(-0.25, 0.25)) << "\n";
std::cout << acos(itv(-0.25, 0.25)) << "\n";
std::cout << atan2(itv(1.), itv(1.)) << "\n";
std::cout << sinh(itv(-0.25, 0.25)) << "\n";
std::cout << cosh(itv(-0.25, 0.25)) << "\n";
std::cout << tanh(itv(-0.25, 0.25)) << "\n";
std::cout << asinh(itv(-0.25, 0.25)) << "\n";
std::cout << acosh(itv(1.5, 2.)) << "\n";
std::cout << atanh(itv(-0.25, 0.25)) << "\n";
// integer power
std::cout << pow(itv(2., 3.), -2) << "\n";
// general power
std::cout << pow(itv(2., 3.), itv(2., 3)) << "\n";
```
The results from the above codes are verified to including the correct results in a strict mathematical sense. However, these mathematical functions return not always the smallest intervals, and therefore, their widths often correspond to a number of floating-point numbers.

### 5.6 Equality and inequality

The following codes works:
```    std::cout << (x < y) << "\n";
std::cout << (x < 1.) << "\n";
std::cout << (1. < x) << "\n";
std::cout << (itv(1.) == 1.) << "\n";
std::cout << (itv(1., 2.) != 3.) << "\n";
std::cout << ("0.20001" > x) << "\n";
```
"true" is returned only when the all points in an interval satisfy a target equality or inequality.

One can get the lower and upper bounds of an interval, and re-write one of them as:
```    // access to endpoints
std::cout << x.lower() << "\n";
std::cout << x.upper() << "\n";
x.lower() = 3.5;
```

### 5.8 division_part1 and division_part2

By default, x / y returns an error if y includes 0. However, when ±∞ is accepted to be bounds of an interval, for example, 1/[-1,1] may return [-∞,-1]∪[1,∞].

"division_part1" and "division_part2" has functions for performing such a division. It is possible that such a division returns two sets of intervals. Therefore, division_part1 returns a bool number; "true" asserts the result from a division is disconnected, i.e., has two sets of intervals. When division_part1 returns "true", the other interval can be computed by division_part2. This procedure is the same as boost.interval. The following is an example code.

```    // division_part1, division_part2
// calculate X / (Y ∖ 0). the result may be divided into two part.
bool parted;
std::cout << division_part1(itv(1., 2.), itv(-3., 4.), parted) << "\n";
// if the result has division_part2, parted is set to true.
if (parted) std::cout << division_part2(itv(1., 2.), itv(-3., 4.)) << "\n";

```
Only when dividing by [0,0], an error is returned.

### 5.9 Other functions

In addition, KV library has the following functions:
```    // static functions
std::cout << itv::whole() << "\n";
std::cout << itv::hull(1., 2.) << "\n";
std::cout << itv::hull(1., z) << "\n";
std::cout << itv::hull(x, y) << "\n";
```
"whole" returns (-∞, ∞).
"hull" returns the convex hull of input arguments.
```    std::cout << width(z) << "\n";
std::cout << median(z) << "\n";
std::cout << mid(z) << "\n";
```
"width" returns the width of an interval, which computed by (upper bound - lower bound) with the up-rounding.
"rad" returns the radios of an interval, which computed by (upper bound - lower bound)/2 with the up-rounding
"median" and "mid" return an (approximate) center of an interval, which computed as (upper bound + lower bound)/2 without controlling a rounding mode. The type of the output is T in the interval.
```    std::cout << norm(z) << "\n";
std::cout << mag(z) << "\n";
std::cout << abs(itv(2., 3.)) << "\n";
std::cout << mig(itv(-2., 1.)) << "\n";
```
"abs" returns the absolute value of an input. If the input is an interval I, abs(I)={abs(x) | x ∈ I}.
"mag" and "norm" both return max{abs(x) | x ∈ I}; the type of the output is T.
"mig" returns min{abs(x) | x ∈ I}; the type of the output is T.
```    std::cout << in(3.9, z) << "\n";
std::cout << subset(x, y) << "\n";
std::cout << proper_subset(y, z) << "\n";
std::cout << overlap(x, y) << "\n";
std::cout << intersect(y, z) << "\n";
```
"in(x, I)" returns "true" if x ∈ I, where the type of x is T and the type of I is interval<T>.
"subset(I, J)" returns "true" if I ⊂ J, where the type of I and J is interval<T>
"proper_subset(I, J)" returns "true" if I is a proper subset of J, where the type of I and J is interval<T>
"overlap(I, J)" returns "true" if I ∩ J ≠ ∅.
"intersect(I, J)" returns I ∩ J, which works only when overlap(I, J) = "true" (KV-library does not supports the empty set).
```    // numeric constants
std::cout << kv::constants<itv>::pi() << "\n";
std::cout << kv::constants<itv>::e() << "\n";
std::cout << kv::constants<itv>::ln2() << "\n";
```
KV-library has the data of constants pi, e, and loge2 a priori.

```    midrad(I, m, r);
```
calculates (approximate) midpoint m and radius r simultaneously such that the ball B(m, r) includes I. If you simply calculate
```    m = mid(I);
```
then B(m,r) may not include I due to rounding error of m.

### 5.11 max, min

```    I3 = max(I1, I2);
I4 = max(I1, I2);
```
"max([a,b],[c,d])" returns [max(a,c),max(b,d)]. "min([a,b],[c,d])" returns [min(a,c),min(b,d)].

## 6. Other interior types (except for double)

This section is incomplete, but the following functions are sufficient for correct interval arithmetics:
• the function for controlling rounding modes (like rdouble.hpp for the type of double)
• the four arithmetic operations
• the sqrt function
• relational operators
• overloading of the ostream (required, e.g, to make "cout << x;" work)
• abs, floor, and frexp
• epsilon, infinity, max, and min in numeric_limits

## 7. Customizing interval arithmetic library

See Customizing interval arithmetic library (in Japanese).

## 8. Emulation of changing rounding modes

When compiling a program, which includes rdouble.hpp, with the option -DKV_NOHWROUND, a computation can be done without any actual changes of rounding modes using the functions of "twosum" and "twoproduct", that is, the all changes of rounding modes are emulated. This improves the portability of this library to various hardwares. For a program including rdd.hpp, the same compile option is also available.

Remark that the above function does NOT work on the 32bit mode on an Intel's CPU, because the function is designed to be able to work on the IEEE754 standard. The 32bit mode on an Intel's CPU does not conform to the IEEE754 standard, and compute floating-point values with the 80bit precision using the FPU.

## 9. Others

• For the type of an interval "::base_type" returns the type of the lower/upper bound of the interval. For example, when "typedef kv::interval<double> itv" is defined, "itv::base_type" returns double.
• kv library does not support the empty set, whereas boost.interval supports it. The author will be grad if one tells him a good way to treat the empty set with preventing a bad influence for a computation performance.