Last update: Sep. 5, 2018

Interval Arithmetic Library

Masahide Kashiwagi

- 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.

- 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.

(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.

#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

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.);

To be compatible with boost.interval, the form likeRemark.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.

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.

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.

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.

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) .

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.

// access to endpoints std::cout << x.lower() << "\n"; std::cout << x.upper() << "\n"; x.lower() = 3.5;

"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.

// 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

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.

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)].

- 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

- 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.