Interval arithmetic is a mathematical technique that handles ranges of values rather than precise numbers. Instead of working with a single value, calculations are performed with intervals representing the range of possible values. Each interval has a lower and upper bound, representing the minimum and maximum possible values respectively.
Contents
You can find examples of interval arithmetic template classes for the C++ programming language but none that are easy to understand and implement. That is why this paper was written “Practical Implementation of Interval Arithmetic”. The paper highlights implementing a general interval template class supporting the float and double type of the C++ programming language.
The primary advantage of interval arithmetic is its ability to provide guaranteed bounds on results, which is especially useful for dealing with rounding errors, measurement uncertainties, or imprecise inputs. When calculations are performed, the results are also intervals that encompass all possible outcomes. This is particularly valuable in fields like numerical analysis, computer science, and engineering, where understanding the precision and reliability of computations is crucial. It allows for more robust error handling and uncertainty management in complex calculations.
Building an interval template software class that can handle all interval arithmetic for IEEE754 floating point is simple math. This paper describes the practical implementation aspect of building an interval template class and discusses the tricks and math behind the implementation. The implementation spans from the basic operators like +,-,*,/ to elementary functions like sqrt(), log(), pow(), exp() and trigonometric functions of sin(), cos(), tan(), asin(), acos() and atan() and hyperebolic functions of sinh(), cosh(), tanh(), asinh(), acosh(), atanh(). Of course, it also supports interval constants like π, ln2, and ln10. The interval arithmetic supports the basic C++ types of float and double.
As usual, we look at the practical implementation aspect of the template class and not so much the theoretical rationale behind it. See [4] Interval Arithmetic: From Principles to Implementation and [5] for some deeper insight into interval arithmetic and the use of IEEE754 floating point format.
This article will be limited to the phase one implementation. The full source and document can be found on my web site at Interval Arithmetic (hvks.com).
- 13th March, 2024. This is a completely rewritten second version of the original form 2015.
The history of interval arithmetic dates back to the mid-20th century and is intertwined with the development of computer science and numerical analysis. Here’s a brief overview of its evolution:
Early Concepts around the 1950s. The roots of interval arithmetic can be traced to the work of mathematicians like Archimedes. However, the formal development began in the 1950s with the advent of digital computers.
The development of interval arithmetic as a formalized field is largely credited to Ramon E. Moore. Moore, while working at the Applied Physics Laboratory at Johns Hopkins University and later at Stanford University, developed the basic rules and operations of interval arithmetic. His work was aimed at addressing the problem of error bounds in calculations, which was crucial for reliable scientific computing.
Expansion and Formalization of interval arithmetic happens in the 1970s-1980s. The field expanded significantly during this period. Researchers began to formalize the arithmetic and explore its theoretical properties. They also started to apply it to various problems in scientific computing, engineering, and error analysis.
With the rise of computer science, interval arithmetic started to be integrated into computer algorithms and software in the 1980s-1990s time frame. This was facilitated by the increasing need for reliable numerical computations in engineering and the sciences.
After the 2000s Interval arithmetic has continued to grow and find new applications. It’s used in fields ranging from numerical analysis to robotics and control systems. The development of specialized software and hardware to efficiently perform interval arithmetic has also been a focus.
Throughout its history, interval arithmetic has been driven by the need for accurate and reliable numerical computations, especially in the face of uncertain data or rounding errors. Its development reflects a broader trend in mathematics and computer science towards more robust and error-resistant methods. As a curious example
Archimedes is one of the earliest known figures to use a form of interval arithmetic in estimating the value of π. He did this by inscribing and circumscribing polygons around a circle and calculating their perimeters. By increasing the number of sides of the polygons, he was able to narrow the bounds and approximate π more closely. At that time, Archimedes couldn't calculate the exact value of π. He determined that π was between (approximately 3.1429) and (approximately 3.1408). This method of approximation was quite sophisticated for its time and is an early historical example of the principles behind interval arithmetic, where values are confined within an upper and lower bound to provide a range of possible values.
Interval arithmetic has several practical applications across various fields due to its ability to handle uncertainties and provide guaranteed bounds on calculations. Here are some of its primary applications starting with Numerical Analysis which ensures accuracy and reliability in computations by accounting for rounding errors and other numerical uncertainties.
In rendering computer graphics and geometric computations, interval arithmetic helps in handling inexact data and calculations, leading to more robust algorithms. For systems that require stability and reliability (control systems), interval arithmetic can manage uncertainties in measurements and model parameters. In planning and control of robotics, it’s used to account for the imprecision in sensor data and to ensure safe movements and decisions. In engineering design, Engineers use interval arithmetic for designing systems with tolerances and variations in material properties. For simulations and modeling in scientific computing where input data might be uncertain or vary over a range, interval arithmetic provides a way to ensure results account for all possible input values. This also extends to optimization problems, error analysis, and propagation which helps in understanding how errors in input data can affect the outcome of calculations, which is crucial in sensitive applications like aerospace and medical equipment.
By providing a framework to handle and propagate uncertainties, interval arithmetic enhances the robustness and reliability of computational tasks in these and other areas.
To bound rounding errors when performing floating point arithmetic, we can use interval arithmetic to keep track of rounding errors. After a series of operations using basic operators like +,-,*, and /, we end with an interval instead of an approximation of the result. The interval width represents the uncertainty of the result but we would know for sure that the correct result will be within this interval.
In interval arithmetic, an interval is typically represented as [a, b], where 'a' is the lower bound and 'b' is the upper bound. Conventionally, it is assumed that a ≤ b, defining what is known as a proper interval. Conversely, if a > b, the interval is termed an improper interval. According to the IEEE 1788 standard for interval arithmetic, both proper and improper intervals are valid and must be correctly handled in computations, as we cannot always assume that a ≤ b. This necessitates determining the lower bound (infimum) of an interval by taking min(a, b), and the upper bound (supremum) by taking max(a, b). This requirement introduces additional steps in the implementation: rather than assuming a ≤ b, code must incorporate functions to find the min and max values for any given interval.
Additionally, an interval where a = b is known as a singleton interval.
For simplicity in formulas and explanations, we often present intervals as if they are proper (a ≤ b). However, in practice, especially in coding implementations, it's crucial to handle intervals correctly according to their proper or improper nature (a>b).
Moreover, it's important to define closed intervals [a, b], semi-closed (or semi-open) intervals (e.g., [a, b) or (a, b]), and open intervals (a, b). Each of these types has different implications for the inclusion or exclusion of the endpoint values in the interval.
The five different interval types are:
| Interval Type | Definition | Symbol | Comments |
| Close | a≤x≤b | [a,b] | |
| Left Open | a<x≤b | (a,b] | Same as right close |
| Right Open | a≤x<b | [a,b) | Same as left close |
| Open | a<x<b | (a,b) | |
| Empty | | ∅ | Empty Interval |
Besides the ∅ (empty interval), the 4 other different types are not strictly needed by the standard. However, it is quite convenient to have them otherwise defined an open or semi-open interval will be challenged. Example, creating an open interval (0,1) will be something like this in C++:
interval<double> i(nextafter(0,+INFINITY), nextafter(1,-INFINITY));
A rather inconvenient and cumbersome way of creating this interval while in my implementation can be done like this:
interval<double> i(0,1,OPEN);
This approach ensures a comprehensive understanding and correct implementation of interval arithmetic, adhering to the IEEE 1788 standard.
An interval is presented by two numbers representing the lower and upper range of the interval:
[a,b] Where ![Image 1]()
(for proper intervals))
The 4 basic operators +,-,*, and / for a proper interval are then defined by:
[a,b] + [c,d] = [a+c,b+d]
[a,b] – [c,d] = [a-d, b-c]
[a,b] * [c,d] = [min(ac,ad,bc,bd),max(ac,ad,bc,bd)]
[a,b] / [c,d] = [min( and [c,d] does not contain 0.
The division can also be rewritten to:
[a,b] / [c,d] = [a,b] * (1/[c,d]) = [a,b] * [1/d,1/c]
Now what we don’t like is the extra work we would need to do for multiplication of intervals. According to the definition, it will require 8 multiplications to do an interval multiplication. We can however get by with a lot less.
Let’s first define an interval class as follows:
| Class if [a,b] | Conditions a≤b |
| Zero | a=0 ^ b=0 |
| Mixed | a<0 ^ b≥0 |
| Positive | a>0 |
| Negative | b<0 |
Instead of the 8 multiplications, we can get by with 2 in most cases and 4 multiplications in worst case based on the following table:
| Class of [a,b] | Class of [c,d] | Multiplication |
| Positive | Positive | [a*c,b*d] |
| Positive | Mixed | [b*c,b*d] |
| Positive | Negative | [b*c,a*d] |
| Mixed | Positive | [a*d,b*d] |
| Mixed | Mixed | [min(a*d,b*c),max(a*c,b*d)] |
| Mixed | Negative | [b*c,a*c] |
| Negative | Positive | [a*d,b*c] |
| Negative | Mixed | [a*d,a*c] |
| Negative | Negative | [b*d,a*c] |
| Zero | - | [0,0] |
| - | Zero | [0,0] |
From a practical standpoint, we can’t just use the above definition “as is”. When performing any basic operations, we will always introduce some rounding errors when using a finite representation like the IEEE754 standard for 32-bit binary floating-point arithmetic and 64-bit binary floating-point arithmetic which is common in most programming languages and CPUs from most vendors like Intel, AMS, etc.
General speaking when performing any of the basic operations like +,-,*,/ we have an error that is bound by the ½β1-t using the default rounding mode in IEEE754 or when rounding towards positive infinity or negative infinity. In IEEE754 binary arithmetic β=2 and t=24 (23+1) for 32-bit float and t=53 (52+1) for 64-bit float, this gives a relative error when performing any of the basic operations
| IEEE754 | Rounding | Towards +∞ or -∞ |
| 32-bits float | ~5.96-8 | ~1.19-7 |
| 64-bits float | ~2.22-16 | ~1.11-16 |
When performing the interval operations, we would need to also control the rounding mode to ensure that the result is really within the interval. 
means rounding towards -∞ and means rounding towards +∞
From an IEEE1788 standard, it is also required that an interval can be empty. It is designated with the symbol ∅. There are a few rules for the empty interval.
[a,b]+∅=[a,b], same goes for ∅+[a,b]:=[a,b]
[a,b]-∅:=[a,b], ∅-[a,b]:=[-b,-a]
[a,b]*∅:=∅, ∅*[a,b]:=∅
[a,b]/∅:=∅, ∅/[a,b]:=∅
Lastly, we also need to be aware of what happens when you square an interval e.g. [a,b]2. If you take the interval [-1,1] and just multiply it with itself using the multiplication formula you get a new interval of [-1,1] which is unnecessarily large. If you look at the plot of x2 in the range [-1,1] you get:
As you can see the result goes from +1 down to 0 and then rises again to +1. The right result for such an operation is [0,1]. The IEEE 1788 standard kind of foreseen this issue and required a compliant implementation to have a function called sqr() for the squaring of intervals. This issue persists for even power of an interval which needs to be handled properly in an implementation.
If you look at the plot for x3, you get:
As you can see, the value of the operation runs through the interval of [-1,1].
The default rounding mode in the IEEE754 floating-point standard is round to nearest (RN for short). However, as shown in the previous section, we need access to round down towards -∞ (RD) and round up towards +∞ (RU).
In interval arithmetic, rounding is a crucial aspect, and it can be implemented either using hardware-based methods or software-based emulation. Hardware-based rounding offers speed and accuracy but lacks flexibility and portability. In contrast, software emulation is more portable and flexible but can be slower and potentially less accurate.
Here are some of the pros of using hardware-based rounding.
Hardware-based rounding is typically faster as it utilizes the processor’s built-in capabilities. It often offers more precise and consistent rounding due to the direct implementation in the CPU. Reduces the overhead of software emulation, making computations more efficient.
However, the advantages in speed come at a cost, since it relies on specific hardware features, which might not be available or consistent across different platforms or processors.
Hardware implementations are harder to customize or update compared to software. Software relying on specific hardware features may face portability issues across different systems.
As required in IEEE754 round to -∞ (RD) and round to +∞ (RU) are available and it is just a simple matter to switch between these modes from the default round to the nearest.
To control the rounding process, the IEEE754 offers at least 4 ways of rounding.
- Rounding to the nearest which is the default rounding method
- Rounding down towards -∞
- Rounding up towards +∞
- Rounding towards zero (truncating)
Depending on the CPU architecture different instructions enabled these different modes. However, it is unfortunately not portable from architecture to architecture. We will just rely on that somehow the following three functions are available that enable these modes.
Fp_near();
Fp_Down();
Fp_up();
As an example of how it can be used. We can implement the interval addition of two intervals [a,b] and [c,d] as follows.
Fp_up();
f=b+d;
Fp_down();
e=a+c;
Fp_near();
return [e,f];
It seems advantageous to use hardware-controlled rounding, however, in [7], they noted that switching back and forth between the different hardware rounding modes is a relatively expensive operation and surprisingly it turns out that using the software rounding method provides a higher performance in many cases. Although it depends highly on the actual CPU architecture.
Software-based rounding for sure has its advantages. Software-based rounding is more portable as it does not depend on specific hardware features. It is easier to customize or adapt to specific needs, such as implementing non-standard rounding modes and it ensures consistent behavior across different platforms and hardware.
The drawback of software-based rounding is that it is typically slower than hardware implementations due to the overhead of emulation and it can be more complex to implement and maintain its accuracy and consistency compared to hardware implementations, depending on the quality of the software emulation.
The software-based rounding takes advantage of round to nearest is the default rounding mode in IEEE754. and two algorithms are useful. The two-sum and two-product algorithms are fundamental in interval arithmetic for maintaining precise bounds. These operations are essential for addressing the problem of rounding errors in floating-point arithmetic. see [7]
- The Two-Sum Algorithm is used to add two floating-point numbers while capturing the error in the result due to rounding. The two-sum operation returns two numbers: the exact sum (as close as the floating-point representation allows) and the error.
- The Two-Product Algorithm is similar to the two-sum, the two-product algorithm multiplies two floating-point numbers, providing the product and the error.
To see how we can take advantage of the two algorithms, let's consider the sum of the arithmetic operation a+b. The result will be rounded to the nearest result (RN). The round to the nearest is either the same as the round down (RD) or the same as rounded up (RU). To figure out which we use the error from the two-sum algorithm. If the error is negative, then:
RN(a+b)==RU(a+b) and RD(a+b) is predecessor(RN(a+b))
If e is positive then:
RN(a+b)==RD(a+b) and RU(a+b) is successor(RN(a+b))
The successor returns the next representable floating-point number towards +∞ and the predecessor returns the next representable floating-point number towards -∞.
The same analogy can be used for multiplication. If the error is negative then:
RN(a*b)==RU(a*b) and RD(a*b) is predecessor(RN(a*b))
If e is positive then:
RN(a*b)==RD(a*b) and RU(a*b) is successor(RN(a*b))
You can extend that to division, reciprocal, and square root, see the table below.
| Operation | Error<0 | Error>0 |
| RN(a+b) | RU=RN, RD=predecessor | RD=RN, RU=successor |
| RN(a-b) | RU=RN, RD=predecessor | RD=RN, RU=successor |
| RN(a*b) | RU=RN, RD=predecessor | RD=RN, RU=successor |
| RN(a/b) | RU=RN, RD=predecessor | RD=RN, RU=successor |
| RN(1/b) | RU=RN, RD=predecessor | RD=RN, RU=successor |
| RN(sqrt(a)) | RU=RN, RD=predecessor | RD=RN, RU=successor |
Two-sum algorithm was invented in 1965 by O. Møller and comes in two variations. The two-sum and the fast two-sum algorithm and implementation of this are below.
The original two-sum algorithm in C++.
std::pair<_IT, _IT> two_sum(_IT a, _IT b)
{
const _IT sum = (a + b);
const _IT a1 = sum - b;
const _IT b1 = sum - a;
const _IT da = a - a1;
const _IT db = b - b1;
const _IT err = da + db;
return std::make_pair(sum, err);
}
This algorithm requires approximately 6 floating-point operations.
A faster version, surprisingly called the fast two-sum algorithm is shown below in C++.
std::pair<_IT, _IT> fasttwo_sum(_IT a, _IT b)
{
const _IT sum = a + b;
_IT err, tmp;
if (abs(a) > abs(b))
{
tmp = sum - a; err = b - tmp;
}
else
{
tmp = sum - b; err = a - tmp;
}
return std::make_pair(sum, err);
}
which only need 3 floating point operations and comparison operations. However, it requires that a≥b for the sum of a+b.
The two products were invented by Rump, see C++ implementation below.
_IT split(_IT a)
{
const _IT C = a*double((1 << 27) + 1);
_IT xhigh, xlow;
xhigh = C - (C - a);
xlow = a - xhigh;
return std::make_pair(xhigh, xlow);
}
std::pair<_IT, _IT> two_prod(_IT a, _IT b)
{
const std::pair<_IT, _IT> x= split(a);
const std::pair<_IT, _IT> y = split(b);
const _IT p = a * b;
const _IT t1 = -p + x.first * y.first;
const _IT t2 = t1 + x.first * y.second;
const _IT t3 = t2 + x.second * y.first;
const _IT err = t3 + x.second * y.second;
if (abs(p) > ldexp(1, -1021)) err = x.second * y.second - ((((p * 0.5 - (x.first * 0.5) * y.first) * 2)
-x.second * y.first) - x.first * y.second);
else
err=x.second*y.second-(((p-x.first*y.first)-x.second*y.first)
-x.first*y.second);
return std::make_pair(p,err);
}
However, as pointed out by [7], if you have access to a fma operation (Fused-Multiply-Add), you can simplify the two-product codes to simple:
The two-product using fma:
std::pair<_IT, _IT> fasttwo_prod(_IT a, _IT b)
{
const _IT p = a * b;
const _IT err = fma(a, b, -p);
return std::make_pair(p, err);
}
Software libraries and systems that support interval arithmetic often implement these algorithms to ensure that all possible values within the intervals are accurately accounted for. Some well-known software and libraries include INTLAB: A MATLAB toolbox for interval arithmetic that provides reliable computations with rounding error bounds. The Boost Interval Arithmetic Library: Part of the Boost C++ Libraries, it supports interval arithmetic operations with various policies for controlling the rounding and handling of interval bounds.
These tools and libraries are used in various fields, from scientific computing to engineering and finance, wherever accurate error tracking and robust interval calculations are needed. They implement not just two-sum and two-products but a whole range of arithmetic and mathematical operations adapted for interval arithmetic.
With the definition of interval arithmetic done and the two crucial algorithms two-sum and two-product, we can now turn to the requirement of an interval class.
It is always important to have your requirements ready before starting the C++ coding process. I usually divide it up into phases so I start with a base implementation and then expand it further through the defined phases. For each phase, you add testing of the phase functionalities.
Phase 1: establishing the base functionality
All the essential operators: +,-,*,/,=,+=,-=,*=,/=.
Constructors for a float, double, or long double to interval type.
Conversion operators from interval to regular float, double, long double, short, int, long.
Support of closed, semi-closed, semi-open, and open intervals.
Support the empty interval set ∅.
cout and cin operator.
Comparison operators ==,!=,<,<=,>,>=.
Essentials methods:
negate: Negate the interval (same as monadic -).
toString: return the string representation of an interval number.
width: return the width of the interval. (standard names this wid).
center: return the center of the interval (standard names this mid).
radius: return the radius of the interval (standard names this rad).
infimum: return the infimum of the interval (standard name it inf).
supremum: return the supremum of the interval. (standard name it sup).
leftinterval: return or set the lower range of the interval as is.
rightinterval: return or set the upper range of the interval as is.
isEmpty: return the Boolean value of emptiness of the interval.
isProper: return true if the interval is a proper interval otherwise false.
isImproper: return true if the interval is an improper interval otherwise false.
isPoint: return true if the interval is a singleton, otherwise false.
isClass: return the classification of the interval.
abs: return the absolute value of the interval.
sqrt: return the square root of the interval.
sqr: return the square of the interval (x2).
Phase 2: expanding the base functionality
mixed data type for arithmetic operations.
manifest interval constants: LN2, LN10, e, PI
exponential (exp) and logarithm functions. (log,log10)
power x<sup>y</sup>
Union: Union of two intervals
Intersection: intersections of two intervals
Interior: is the interval an interior to another interval
Precedes: Does the interval precede another interval
Subset: is the interval a subset of another interval
Inclusion: is the interval included in another interval
In: is the float or double number within the interval
Phase 3: Completed all the auxiliary functions
Trigonometric interval functions
Hyperbolic interval functions:
Phase 4:
Expand support for arbitrary precision interval type
other useful functions e.g. gamma, beta, etc.
Is all about establishing the base functionality of the interval class.
With the implementation split into phases, we are now ready to begin our initial design.
We need to create a template class for our interval functionality. In essence, similar to the standard complex template library in C++. Also, the class needs to hold the left and right sides of the interval plus the interval type to be able to cope with Closed, Semi-open, or fully open intervals as previously discussed. Not surprisingly, we end up with this template class definition and the private parts and variables.
template<class _IT> class interval {
_IT left; _IT right; enum interval_type type; RIGHT_OPEN, (RIGHT_CLOSE is synonym for
LEFT_OPEN and LEFT_CLOSE same as RIGHT_OPEN
public:
};
We have also established our two functions, fasttwo_sum() and fasttwo_prod() as internal private functions to the class and we get the following:
template<class _IT> class interval {
_IT left; _IT right; enum interval_type type;
_IT split(_IT a)
{
const _IT C = a * double((1 << 27) + 1);
_IT xright, xleft;
xright = C - (C - a);
xleft = a - xright;
return std::make_pair(xright, xleft);
}
static std::pair<_IT, _IT> fasttwo_sum(_IT a, _IT b)
{
const _IT sum = a + b;
_IT err, tmp;
if (abs(a) > abs(b))
{
tmp = sum - a; err = b - tmp;
}
else
{
tmp = sum - b; err = a - tmp;
}
return std::make_pair(sum, err);
}
std::pair<_IT, _IT> fasttwo_prod(_IT a, _IT b)
{
const _IT p = a * b;
const _IT err = fma(a, b, -p);
return std::make_pair(p, err);
}
public:
};
That was the private part of our template class. Now to the constructors for the class.
As our requirement goes, our interval class can be declared with zero, one, two, or three parameters. Like in:
Interval<double> i;
Interval<float> fi(2.0)
Interval<double> di(1.0,2.0);
Interval<double> dit(3.0,4.0,OPEN);
Interval<float> dit2(3,2);
interval(); interval(const _IT& d); interval(const _IT& l, const _IT& h, const enum interval_type t=CLOSE);
The interval type is defined below:
enum interval_type { EMPTY, CLOSE, LEFT_OPEN, RIGHT_OPEN, OPEN,
RIGHT_CLOSE=LEFT_OPEN, LEFT_CLOSE=RIGHT_OPEN };
We would also need some coordinating functions. Example, get the left or right value of the interval or change them for that matter. Likely we need to find the width (wid) of the interval, the radius (rad) of the interval, and the midpoint (mid) for converting an interval to a scalar element and we also need to define the required infimum (inf) and supremum (sup) plus the “is” method functions as required by the IEEE 1788 standard.
We also need some direct method to get and set the left and right interval plus getting or setting the interval type ([], (], [), []) and I use leftinterval() and rightinterval() for the interval ends:
_IT rightinterval() const; _IT leftinterval() const; _IT rightinterval(const _IT&); _IT leftinterval(const _IT&); enum interval_type intervaltype() const; enum interval_type intervaltype(enum interval_type);
_IT inf() const; _IT sup() const; _IT mid() const; _IT rad() const; _IT wid() const; _IT mig() const; _IT mag() const;
bool isEmpty() const;
bool isPoint() const;
bool isImproper() const;
bool isProper() const;
Previously, I introduced the term interval class as a way to classify a given interval type. Example, zero, positive, negative, mixed, etc. This can also be a helpful function to support. And the toString() function that generates a string representation of the interval.
enum int_class isClass() const;
std::string toString() const;
the enumerated type int_class is defined as follows:
enum int_class { NO_CLASS, ZERO, POSITIVE0, POSITIVE1, POSITIVE, NEGATIVE0,
NEGATIVE1, NEGATIVE, MIXED };
Is all about defining conversion operators to convert an interval to any of the built-in types in C++.
operator short() const operator int() const operator long() const operator long long() const operator unsigned short() const operator unsigned int() const operator unsigned long() const operator unsigned long long() const operator double() const; operator float() const;
This is all the C++ assignments operators, see below.
interval<_IT>& operator= ( const interval<_IT>& );
interval<_IT>& operator+=( const interval<_IT>& );
interval<_IT>& operator-=( const interval<_IT>& );
interval<_IT>& operator*=( const interval<_IT>& );
interval<_IT>& operator/=( const interval<_IT>& );
These are the +,-,*,/ and the monadic version of +,- plus the Boolean ==,!=,>,<,>=,<= operators.
template<class _IT> interval<_IT> operator+( const interval<_IT>&, const interval<_IT>& );
template<class _IT> interval<_IT> operator+( const interval<_IT>& ); template<class _IT> interval<_IT> operator-( const interval<_IT>&, const interval<_IT>& );
template<class _IT> interval<_IT> operator-( const interval<_IT>& ); template<class _IT> interval<_IT> operator*(const interval<_IT>&, const interval<_IT>&);
template<class _IT> interval<_IT> operator/( const interval<_IT>&, const interval<_IT>& );
with the class declaration defined and the outer declaration in place, we can now move to the implementation of it.
Most of the implementation is simple and follows the class definition.
- Constructors
- Coordinate methods
- Conversion operators
- Essential operators
The constructors are simple and just populate the interval class with initial values. If no parameters are given, the interval will be initialized as an EMPTY (∅) interval. Notice that for the three arguments constructor, the third argument is defaulted to a closed interval, if not present.
template<class _IT> inline interval<_IT>::interval()
{ left = _IT(0); right = _IT(0); type = EMPTY; }
template<class _IT> inline interval<_IT>::interval(const _IT& d)
{ left = _IT(d); right = _IT(d); type = CLOSE; }
template<class _IT> inline interval<_IT>::interval
(const _IT& l, const _IT& h, const enum interval_type t)
{ left = l; right = h; type = t; }
These are all the methods as defined in the interval class. Most of these methods are relatively easy to implement. For the methods, we have:
rightinterval()leftinterval()intervaltype()
These are simple functions that just return the left, right or type of interval. There is no computation needed and it just returns the current values. There is also an overloaded version that allows you to set these private variables. There is also the isClass() and toString() methods.
The next group is the required IEEE1788 standard functions:
inf()sup()mid()rad()wid()mig()mag()isEmpty()isPoint()isImproper()isProper()
Two methods inf() and sup() requires a little explanation.
The inf() function stands for "infimum" or "greatest lower bound (GLB)." It returns the lowest possible value in an interval, effectively representing the lower bound of the interval. This function is crucial for ensuring that any operations or comparisons involving intervals take into account the lowest possible value that an element of the interval can assume, thus maintaining accuracy and reliability in mathematical computations. This mean that regardless of whether the interval is proper (left≤right) or improper (left>right), it still returns the infimum of the interval. This is also valid regardless of whether the interval is closed, open, or semi-open.
The sup() function stands for "supremum". It returns the highest possible value in an interval, effectively representing the upper bound of the interval. Similarly, the sup() function ensures that the highest possible value within an interval is considered during operations or comparisons. This is vital for encompassing the full range of possible values within an interval and for accurate interval arithmetic computations. Again, regardless of whether the interval is proper (left≤right) or improper (left>right), it still returns the infimum of the interval. This is also valid regardless of whether the interval is closed, open, or semi-open
In interval arithmetic, these functions are used to define intervals (e.g., [inf(x), sup(x)] for an interval x) and to perform arithmetic operations between intervals. When operations such as addition, subtraction, multiplication, or division are performed between intervals, the inf() and sup() functions help determine the resulting interval by calculating the new bounds based on the operation performed and in consideration of the interval type and properness of the interval.
template<class _IT> inline _IT interval<_IT>::rightinterval() const
{ return right; }
template<class _IT> inline _IT interval<_IT>::leftinterval() const
{ return left; }
template<class _IT> inline _IT interval<_IT>::rightinterval(const _IT& r)
{
if (type == EMPTY)
type = CLOSE;
right = r;
return right;
}
template<class _IT> inline _IT interval<_IT>::leftinterval(const _IT& l)
{
if (type == EMPTY)
type = CLOSE;
left = l;
return left;
}
template<class _IT> inline enum interval_type interval<_IT>::intervaltype() const
{ return type; }
template<class _IT> inline enum interval_type interval<_IT>::intervaltype
(const enum interval_type to)
{
interval<_IT> x = *this;
const _IT infi(INFINITY);
if (x.type != to)
{
if (this->isImproper())
{ std::swap(x.left, x.right);
if (x.type == LEFT_OPEN)
x.type = RIGHT_OPEN;
else
if (x.type == RIGHT_OPEN)
x.type = LEFT_OPEN;
}
switch (x.type)
{
case CLOSE:
if (to == LEFT_OPEN || to == OPEN)
x.left = nextafter(x.left, -infi);
if (to == RIGHT_OPEN || to == OPEN)
x.right = nextafter(x.right, +infi);
break;
case OPEN:
if (to == RIGHT_OPEN || to == CLOSE)
x.left = nextafter(x.left, +infi);
if (to == LEFT_OPEN || to == CLOSE)
x.right = nextafter(x.right, -infi);
break;
case LEFT_OPEN:
if (to == RIGHT_OPEN || to == CLOSE)
x.left = nextafter(x.left, +infi);
if (to == RIGHT_OPEN || to == OPEN)
x.right = nextafter(x.right, +infi);
break;
case RIGHT_OPEN:
if (to == LEFT_OPEN || to == OPEN)
x.left = nextafter(x.left, -infi);
if (to == LEFT_OPEN || to == CLOSE)
x.right = nextafter(x.right, -infi);
break;
}
x.type = to;
if (this->isImproper())
{
std::swap(x.left, x.right);
}
*this = x;
}
return x.type;
}
template<class _IT> inline _IT interval<_IT>::inf() const
{
if (type == CLOSE)
return min(left, right);
_IT adjustedLeft = left;
_IT adjustedRight = right;
const _IT infi(INFINITY);
if (type == LEFT_OPEN || type == OPEN)
adjustedLeft = nextafter(left, (left <= right) ? +infi : -infi);
if (type == RIGHT_OPEN || type == OPEN)
adjustedRight = nextafter(right, (left <= right) ? -infi : +infi);
return std::min(adjustedLeft, adjustedRight);
}
template<class _IT> inline _IT interval<_IT>::sup() const
{
if (type == CLOSE)
return max(left, right);
_IT adjustedLeft = left;
_IT adjustedRight = right;
const _IT infi(INFINITY);
if (type == LEFT_OPEN || type == OPEN)
adjustedLeft = nextafter(left, (left <= right) ? +infi : -infi);
if (type == RIGHT_OPEN || type == OPEN)
adjustedRight = nextafter(right, (left <= right) ? -infi : +infi);
return std::max(adjustedLeft, adjustedRight);
}
template<class _IT> inline _IT interval<_IT>::mid() const
{ if (right == left) return left; else return (right + left) / _IT(2); }
template<class _IT> inline _IT interval<_IT>::rad() const
{ _IT r; r = (right - left) / _IT(2); return r; }
template<class _IT> inline _IT interval<_IT>::wid() const
{ _IT r; r = right - left; if (r < _IT(0)) r = -r; return r; }
template<class _IT> inline _IT interval<_IT>::mig() const
{
_IT l = inf();
_IT r = sup();
l = abs(l);
r = abs(r);
return std::min(l, r);
}
template<class _IT> inline _IT interval<_IT>::mag() const
{
_IT l = inf();
_IT r = sup();
l = abs(l);
r = abs(r);
return std::max(l, r);
}
template<class _IT> inline bool interval<_IT>::isProper() const
{ return left<=right; }
template<class _IT> inline bool interval<_IT>::isImproper() const
{ return right<left; }
template<class _IT> inline bool interval<_IT>::isPoint() const
{ return left==right; }
template<class _IT> inline bool interval<_IT>::isEmpty() const
{ return type==EMPTY; }
template<class _IT> inline enum int_class interval<_IT>::isClass() const
{
if (left == _IT(0) && right == _IT(0)) return ZERO;
if (left == _IT(0) && right > _IT(0)) return POSITIVE0;
if (left > _IT(0) && right > _IT(0)) return POSITIVE1;
if (left >= _IT(0) && right > _IT(0)) return POSITIVE;
if (left < _IT(0) && right == _IT(0)) return NEGATIVE0;
if (left < _IT(0) && right < _IT(0)) return NEGATIVE1;
if (left < _IT(0) && right <= _IT(0)) return NEGATIVE;
if (left < _IT(0) && right > _IT(0)) return MIXED;
return NO_CLASS;
}
template<class _IT> inline std::string interval<_IT>::toString() const
{
std::string s;
enum interval_type t = intervaltype();
std::ostringstream strs;
strs.precision(25);
strs << (t == LEFT_OPEN || t == OPEN ? "(" : "[");
strs << left;
strs << ",";
strs << right;
strs << (t == RIGHT_OPEN || t == OPEN ? ")" : "]");
return strs.str();
}
It is trivial to implement the conversion operators since it just has to return the midpoint of the interval when conversion to any integers type or float and double types. If the interval is empty, the conversion is undefined.
template<class _IT> inline interval<_IT>::operator short() const
{ return static_cast<short>(mid());
}
template<class _IT> inline interval<_IT>::operator int() const
{ return static_cast<int>(mid());
}
template<class _IT> inline interval<_IT>::operator long() const
{ return static_cast<long>(mid());
}
template<class _IT> inline interval<_IT>::operator long long() const
{ return static_cast<long long>(mid());
}
template<class _IT> inline interval<_IT>::operator unsigned short() const
{ return static_cast<unsigned short>(mid());
}
template<class _IT> inline interval<_IT>::operator unsigned int() const
{ return static_cast<unsigned int>(mid());
}
template<class _IT> inline interval<_IT>::operator unsigned long() const
{ return static_cast<unsigned long>(mid());
}
template<class _IT> inline interval<_IT>::operator unsigned long long() const
{ return static_cast<unsigned long long>(mid());
}
template<class _IT> inline interval<_IT>::operator double() const
{ return static_cast<double>(mid());
}
template<class _IT> inline interval<_IT>::operator float() const
{ return static_cast<float>(mid());
}
Now that we can control rounding via software emulation, it is pretty straightforward to implement interval arithmetic in a portable way. An assignment is just a copy from the right to the left side of the = operator.
template<class _IT> inline interval<_IT>& interval<_IT>::operator=( const interval<_IT>& a )
{
left = a.left;
right = a.right;
type = a.type;
return *this;
}
The Interval addition is:
Code for addition using the fastwo_sum. Where .inf() returns the lowest number of the interval. And .sup() returns the highest number in the interval.
template<class _IT> inline interval<_IT>& interval<_IT>::operator+=(const interval<_IT>& a)
{
if (a.type == EMPTY)
return *this;
if (type == EMPTY)
return (*this = a);
const _IT infi(INFINITY);
std::pair<_IT, _IT> xlow = fasttwo_sum(inf(), a.inf());
std::pair<_IT, _IT> xright = fasttwo_sum(sup(), a.sup());
left = xlow.first;
right = xright.first;
if (xlow.second < _IT(0))
left = nextafter(left, -infi);
if (xright.second > _IT(0))
right = nextafter(right, +infi);
type = CLOSE;
return *this;
}
The interval subtraction is:
Code for subtraction using the fasttwo_sum.
template<class _IT> inline interval<_IT>& interval<_IT>::operator-=(const interval<_IT>& a)
{
if (a.type == EMPTY)
return *this;
if (type == EMPTY)
return (*this = -a);
const _IT infi(INFINITY);
std::pair<_IT, _IT> xlow = fasttwo_sum(inf(), -a.sup());
std::pair<_IT, _IT> xright = fasttwo_sum(sup(), -a.inf());
left = xlow.first;
right = xright.first;
if (xlow.second < _IT(0))
left = nextafter(xlow.first, -infi);
if (xright.second > _IT(0))
right = nextafter(xright.first, +infi);
type = CLOSE;
return *this;
}
The Interval multiplication is:
Code for multiplication using the fasttwo_product algorithm and with some optimization.
template<class _IT> inline interval<_IT>& interval<_IT>::operator*=(const interval<_IT>& a)
{
if (type == EMPTY)
return *this;
if (a.type == EMPTY)
return (*this = a);
_IT al = inf();
_IT ah = sup();
_IT bl = a.inf();
_IT bh = a.sup();
auto multiply = [&](const _IT& x, const _IT& y)
{
std::pair<_IT, _IT> tmp = interval<_IT>::fasttwo_prod(x, y);
interval<_IT> res;
const _IT infi(INFINITY);
res.left = res.right = tmp.first;
if (tmp.second < _IT(0))
res.left = nextafter(tmp.first, -infi);
if (tmp.second > _IT(0))
res.right = nextafter(tmp.first, +infi);
return res;
};
interval<_IT> itmp;
type = CLOSE;
if (al >= _IT(0) && bl >= _IT(0))
{ itmp = multiply(al, bl);
left = itmp.left;
itmp = multiply(ah, bh);
right = itmp.right;
return *this;
}
if (ah < _IT(0) && bh < _IT(0))
{ itmp = multiply(al, bl);
right = itmp.right;
itmp = multiply(ah, bh);
left = itmp.left;
return *this;
}
if (al >= _IT(0) && bh < _IT(0))
{ itmp = multiply(ah, bl);
left = itmp.left;
itmp = multiply(al, bh);
right = itmp.right;
return *this;
}
if (ah < _IT(0) && bl >= _IT(0))
{ itmp = multiply(al, bh);
left = itmp.left;
itmp = multiply(ah, bl);
right = itmp.right;
return *this;
}
itmp = multiply(al, bl);
left = itmp.left;
right = itmp.right;
itmp = multiply(al, bh);
left = min(left, itmp.left);
right = max(right, itmp.right);
itmp = multiply(ah, bl);
left = min(left, itmp.left);
right = max(right, itmp.right);
itmp = multiply(ah, bh);
left = min(left, itmp.left);
right = max(right, itmp.right);
return *this;
}
The Interval division is:
template<class _IT> inline interval<_IT>& interval<_IT>::operator/=(const interval<_IT>& b)
{
const _IT infi(INFINITY);
if (type == EMPTY)
return *this;
if (b.type == EMPTY)
return (*this = b);
interval<_IT> tmp = *this; auto inverse = [&](const _IT& y, const bool up)
{
_IT res = _IT(1) / y, r;
r = -fma(res, y, _IT(-1));
if (up == false)
{
if (r < _IT(0))
res = nextafter(res, -infi);
}
else
{
if (r > _IT(0))
res = nextafter(res, +infi);
}
return res;
};
_IT bl = b.inf();
_IT bh = b.sup();
if (bl == _IT(0) && bh == _IT(0))
{
left = -infi;
right = +infi;
*this *= tmp;
return *this;
}
if (bh == _IT(0))
{ right = inverse(bl, true);
left = -infi;
*this *= tmp;
return *this;
}
if (bl == _IT(0))
{ left = inverse(bh, false);
right = +infi;
*this *= tmp;
return *this;
}
left = inverse(bh, false);
right = inverse(bl, true);
*this *= tmp;
return *this;
}
We simply compute the dyadic operators with the use of the essential operators to simplify the code.
template<class _IT> inline interval<_IT> operator+(const interval<_IT>& a, const interval<_IT>& b)
{
interval<_IT> result(a);
result += b;
return result;
}
template<class _IT> inline interval<_IT> operator+(const interval<_IT>& a)
{
return a;
}
template<class _IT> inline interval<_IT> operator-(const interval<_IT>& a, const interval<_IT>& b)
{
interval<_IT> result(a);
result -= b;
return result;
}
template<class _IT> inline interval<_IT> operator-(const interval<_IT>& a)
{
interval<_IT> result(0);
result -= a;
return result;
}
template<class _IT> inline interval<_IT> operator*(const interval<_IT>& a, const interval<_IT>& b)
{
interval<_IT> result(a);
result *= b;
return result;
}
template<class _IT> inline interval<_IT> operator/(const interval<_IT>& a, const interval<_IT>& b)
{
interval<_IT> result(a);
if (result == b)
{
enum int_class intclass = b.isClass();
if (intclass != ZERO && intclass != POSITIVE0 && intclass != NEGATIVE0)
{
result = interval<_IT>(1, 1);
return result;
}
}
result /= b;
return result;
}
In this section, we handle all the 6 comparison operators. ==,!=, >=, <=, >, <. Notice that we also need to be able to handle the empty interval ∅. If both intervals are ∅, then we return false. If one of the intervals is ∅, then we return true. Notice that some implementations return undefined when both intervals are ∅.
template<class _IT> inline bool operator==(const interval<_IT>& a, const interval<_IT>& b)
{
if (a.intervaltype() == EMPTY && b.intervaltype() == EMPTY)
return false; if (a.intervaltype() == EMPTY || b.intervaltype() == EMPTY)
return true; about Improper interval
return a.inf() == b.inf() && a.sup() == b.sup();
}
template<class _IT> inline bool operator!=(const interval<_IT>& a, const interval<_IT>& b)
{
return !(a == b);
}
template<class _IT> inline bool operator>=(const interval<_IT>& a, const interval<_IT>& b)
{
if (a.intervaltype() == EMPTY && b.intervaltype() == EMPTY)
return false;
if (a.intervaltype() == EMPTY || b.intervaltype() == EMPTY)
return true;
if (a.inf() >= b.sup())
return true;
return false;
}
template<class _IT> inline bool operator<=(const interval<_IT>& a, const interval<_IT>& b)
{
if (a.intervaltype() == EMPTY && b.intervaltype() == EMPTY)
return false;
if (a.intervaltype() == EMPTY || b.intervaltype() == EMPTY)
return true;
if (a.sup() <= b.inf())
return true;
return false;
}
template<class _IT> inline bool operator>(const interval<_IT>& a, const interval<_IT>& b)
{
if (a.intervaltype() == EMPTY && b.intervaltype() == EMPTY)
return false;
if (a.intervaltype() == EMPTY || b.intervaltype() == EMPTY)
return true;
if (a.inf() > b.sup() )
return true;
return false;
}
template<class _IT> inline bool operator<(const interval<_IT>& a, const interval<_IT>& b)
{
if (a.intervaltype() == EMPTY && b.intervaltype() == EMPTY)
return false;
if (a.intervaltype() == EMPTY || b.intervaltype() == EMPTY)
return true;
if (a.sup() < b.inf())
return true;
return false;
}
We need to add a special code for handling x2. See previous sections. For the sqrt(x), we again can use the fma() library function to quickly bound the interval.
template<class _IT> inline interval<_IT> sqr(const interval<_IT>& x)
{
_IT left = x.inf();
_IT right = x.sup();
_IT tmpl = left;
_IT tmpr = right;
interval<_IT> r(0);
left *= left;
right *= right;
r.rightinterval(max(left, right));
if ( tmpl > 0 && tmpr > 0)
r.leftinterval(min(left, right));
return r;
}
template<class _IT> inline interval<_IT> sqrt(const interval<_IT>& x)
{
const _IT infi(INFINITY);
_IT sq, r;
_IT left, right;
sq = sqrt(x.inf());
left = sq;
r = -fma(sq, sq, -x);
if (r < _IT(0))
left = nextafter(left, -infi);
sq = sqrt(x.sup());
right = sq;
r = -fma(sq, sq, -x);
if (r > _IT(0))
right = nextafter(right, +infi);
return interval<_IT>(left, right);
}
In interval arithmetic, the absolute value (abs) of an interval is defined in a way that captures all possible absolute values of any number within the interval. The result is always a non-negative interval. The definition depends on the position of the interval concerning zero:
If the interval is entirely non-negative (i.e., both its lower and upper bounds are greater than or equal to zero), then the absolute value of the interval is the interval itself, since all values within it are already non-negative.
If the interval is entirely negative (i.e., both its lower and upper bounds are less than zero), then the absolute value of the interval is obtained by taking the absolute values of its bounds and swapping them (since the lower bound will be the most negative and, thus, have the highest absolute value).
If the interval spans zero (i.e., its lower bound is negative and its upper bound is positive), then the absolute value of the interval is from zero to the maximum of the absolute values of the lower and upper bounds.
Source code for abs(x):
template<class _IT> inline interval<_IT> abs( const interval<_IT>& a )
{
if (a.inf() >= _IT(0) ) return a;
else
if (a.sup() < _IT(0) ) return -a;
return interval<_IT>(_IT(0), max(-a.inf(),a.sup()));
}
The standard output cout and input cin are straightforward. cout outputs the interval in the format of:
[ left, right ]
When the interval is closed:
( left, right ] or [ left, right ) for semi-open intervals and
And:
( left, right ) for fully open intervals.
If the interval is empty, an ∅ is outputted.
While cin can handle a little variation to the above. Any syntax for the input is legal.
[ left, right ]
[ interval ] which is equivalent to [ interval, interval ]
you can mix in ‘(‘ or ‘)’ instead of square bracket ‘[‘ or ‘]’ to indicate semi or full-open intervals.
This was the goal for the first phase which defined and implemented all the base base-level functionality that is needed for the next phases.
However, we do have a shortcoming. Where the result is not completely accurate. This has to do with the constructor. The constructor parameter is of the template type _IT (either float, double, or long double). As per C++ standard rules, when an actual parameter is of a lesser type than _IT, then an automatic conversion is applied to the desired _IT type. However, the conversion happens before executing the constructor and therefore any inaccuracies in the conversion are not caught by the interval class. An example will clarify the issue.
Interval<float> x(16777216);
Interval<float> y(16777217);
Will create the same interval [16777216.0,16777216.0] for both x and y. The issue is that the integer 16777217 is larger than what accurately can be converted to a float 23-bit mantissa while 16777216 just fits in.
The same is true for converting a 64-bit double to a 32-bit float type. Even for an interval<double> variable, a 64-bit integer larger than 9,007,199,254,740,992 will not convert correctly to a double (a double mantissa is 52-bit.
In phase 2, we will introduce mixed data type which magically allows us to address and fix the above issue.
This concludes Phase One of the practical implementation of an interval template library. We now have the basic functionality established and we can continue implementing Phase Two and Phase Three.
To read the entire paper which includes the phase two and phase three implementation, the full source and document can be found on my web site at Interval Arithmetic (hvks.com).
- Intel 64 and IA32 Architectures Software Developers Manual. Volume 2. June 2014
- Wilkinson, J H, Rounding errors in Algebraic Processes, Prentice-Hall Inc, Englewood Cliffs, NJ 1963
- Richard Brent & Paul Zimmermann, Modern Computer Arithmetic, Version 0.5.9 17 October 2010; http://maths-people.anu.edu.au/~brent/pd/mca-cup-0.5.9.pdf
- T. Hickey & M.H. van Emden. Interval Arithmetic: from Principles to Implementation
- Hend Dawood, Theories of Interval Arithmetic, Mathematical Foundations and Applications, Lambert Academic Publishing. 2011. ISBN: 978-3-8465-0152-2
- Wikipedia, Interval Arithmetic, https://en.wikipedia.org/wiki/Interval_arithmetic,
- Daniel Pfeffer, Interval arithmetic using round-to-nearest, https://www.codeproject.com/Articles/1040839/Interval-arithmetic-using-round-to-nearest-mode-pa
- N.T. Hayes, February 2013, Standard for Modal Interval Arithmetic.
- J.G. Role, Interval Arithmetic and Analysis: An Introduction
