Final update: March 29, 2016

Interval Arithmetic Library

Masahide Kashiwagi

1. Introduction

2. Features

3. Supported environment

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.

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 loge2 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:

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