Last update: Jan. 6, 2022

Interval Arithmetic Library

Masahide Kashiwagi

1. Introduction

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.hpp (to be introduced later), we need to change the rounding mode. Recent compilers are C99 compilant, and we can change the rounding mode with fesetround in fenv.h provided by C99. In addition, it may be possible to speed up with the compile option, see 20. Changing Rounding Mode and Compile Options 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");

    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 narrowest intervals, and not so fast.

For details of implementation of mathematical functions, see Simple implementation of interval mathemarical functions (Version: 2018/9/5, written in Japanese) .

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.

5.7 Access to bounds

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 << rad(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 log_e2 a priori.

5.10 midrad

    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);
    r = rad(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. Conversion between intervals with different base type

See Conversion between intervals with different base type (in Japanese).

9. Others

It is considered to be able to handle intervals where one or both of the end points are infinite. The implementation details are described in Summary of interval arithmetic with infinity at the endpoint (in Japanese).
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.