Boost.Real a C++ approach to range arithmetic

Documentation

• Project documentation main page
• Doxygen documentation
• Medium post

Information

Source code: https://github.com/BoostGSoC18/Real

Final release of the library during the GSoC 18: https://github.com/BoostGSoC18/Real/releases/tag/v1.0.0-beta
The final release tag serves as a pointer to the final commit made in the Google Summer of Code 2018.

Commits authors: All commits in the GitHub project were made by me within the GSoC 18 program with the exception of the following commits:

Hash: 92b8cd919650c3142115d53bc366ed4d184198d6
Hash: f5618648292a897f23aa607b612f4dfe9d84c752
Hash: 15d76b303fffa1450011a438cadc43a1ee6087ed
Author: Damian Vicino (My mentor in this project)
Purpose: Adding configuration files for Travis CI and codecov support.
Hash: 00ae427f205b9760696178765a24f928ba16d59a
Author: Damian Vicino (My mentor in this project)
Purpose: Set CMakeFile for header only library and testing. 

Project completion:
During the GSoC 18, I have completed the exploratory work predefined for the project and we got a functional library. But even though we have a functional library that can be downloaded and included in other projects, we still need to work on some required performance improvements before submitting the code for inclusion to the Boost C++ library. My mentor and I are looking forward to book a review for the library soon hoping that it will be accepted for the next release of the Boost C++ library.

Introduction

In this blog, I will try to explain the project I have been working for the Google Summer of Code 2018 (GSoC 18). I have applied for the organization Boost C++ to implement a new exploratory library to add to Boost called Boost.Real. The idea is to explore a new abstract representation for real numbers that dynamically set precision by usage context.

The problem addressed by boost::real
Several times, when dealing with complex mathematical calculus, numerical errors can be carried from one operation to the next and after several steps, those errors may significantly increase obtaining an incorrect result. Commonly, numbers representation are limited to a fixed amount of memory (as for example floating point), in these representations, the error is generated because numbers requiring more memory than the available are truncated. Another problem is that error is not tracked, so we don’t even know the magnitude of the final error generated in a complex operation. In C++, the most common representation of reals is by native floating point numbers (Figure 1). These numbers represent finite subsets of the real numbers (which is an infinite set).

Figure 1: Floating point number representation in base 10, a number is represented by a vector of digits (numbers between 0 and 9) and an integer exponent.

The major problem with representing the numbers using a fixed finite memory (as the floating point numbers) is that it is not possible to represent an infinite set of numbers as the set of the real numbers. Instead, the real number line must be discretized to represent only some of them. Because the representation uses an exponent to handle the dot position that divides the integer from the fractional part of the numbers, there are more numbers represented closer to the zero than far from the zero. This discretization is shown in Figure 2.

Figure 2: The dispersion of the represented numbers over the real number line

When an operation results in a number that cannot be represented, an approximation function is used to decide which one of the two closest representable numbers is returned as the operation result. This approximation introduces error to the calculation as shown in Figure 3.

Figure 3: Approximation error generated by the half-precision floating-point representation number representation when the result is not a fully represented number as shown in section 3.1.1 of [5].

Another major problem when dealing with real numbers is the representation of the irrational numbers as for example π or e. These numbers are not handled by native number data types causing limitations when calculations are based on those numbers. One way to deal with this problem is to define a truncation that is good enough for the calculation purposes, but in some cases the error depends in a chain of operations that lead to unbound error or that are unknown in advance to do proper analysis

The solution proposed in Boost.Real
Boost::real is a real number representation data type we have started in the Google Summer of Code 2018 to addresses the mentioned issues using range arithmetic [1]. Range arithmetic is an abstract number representation that improves the conventional floating point representation by adding two new component:

1. It defines the precision as a dynamical attribute to be determined in run-time.
2. It allows error monitoring.

The main goal of range arithmetic is to represent a real number x as a program that returns a finite or infinite set of intervals containing the number x. Then, instead of comparing two numbers using the traditional mathematical operator, we use the operator defined by interval arithmetic to determine which numbers interval is lower than the other.

We define the precision of an interval presenting the number x as the maximum length between any number of the interval and the number x. Thus, given two intervals, the smaller interval has more precision than the biggest interval.

A number x is represented by a program x(k) that return intervals containing the number x. Those intervals are sorted from lowest to higher precision as shown in Figure 4. It is important that the intervals returned by the program x(k) have not the same precision.

Figure 4: x(k) is an interval returned by the program that represent the number x (contained within the interval). The higher it is k, the smaller is the interval and more precision is be obtained by increasing k.

For this purpose, boost::real has a const_precision_iterator that iterates the set of intervals representing the number. Because all the intervals in the set have different precisions and they are sorted from lower to higher precision, each time the iterator is increased (using the operator ++), the number precision increases as well. If in one iteration, the interval is the represented number itself (if possible) then there are no more intervals to iterate. This simple idea is the key concepts of this project and we suggest to take the needed time to correctly understand it before continuing.

Also, to allow representing irrational numbers as π or e, boost::real has a constructor that takes as its parameter a function pointer or a lambda expression that for any given unsigned integer n ≥ 0, the function returns the n-th digit of the represented number. For example, the number 1/3 can be easily represented by a program that for any input n ≥ 0, the function returns 3 and by setting the number exponent as 0.

[] (unsigned int n) -> int { return 3; }

Representing numbers as abstract entities that can be approximated as much as possible allows us to work with them even if they hay an infinite representation that we can’t represent. Using range arithmetic we are now able to use the necessary precision required by the operation, and the precision can be different for each operation. For example, if we want to know if 1/3 is less than 4, we just need to iterate the intervals of the number 1/3 until the upper boundary of the approximation interval is less than 4, but if we need to know if 1/3 is less than 0.33334, then we need to iterate the intervals of the number 1/3 until it is less than 0.33334.

Also, we are able to compare 1/3 with irrational numbers as for example the Champernowne binary number which is composed by successively concatenate the binary numbers as its digits. The Champernowne binary number starts as 0.11011100101110111. Because for the 1/3 number we can have the approximation interval [0.3, 0.4] and for the Champernowne number we can have the approximation interval [0.1, 0.2] then we can determine that Champernowne binary is lower than 1/3 even if both numbers have an infinite representation. This example was obvious but shows the possibility to work with irrational numbers.

In the section following sections (GSoC stage 1, GSoC stage 2 and GSoC stage3) , we will describe the work that we have done during the Google Summer of Code 2018 divided into 3 stages (one for each evaluation). Our intentions in this blog are:

1. Explain the designed architecture of the new data type that we have created.
2. Walk through the design and work evolution during the Google Summer of Code. Because this is a new project started from scratch, a considerable part of the jobs relied on designing the correct strong and clean structures that can be easy to understand, use and contribute by other people.

GSoC stage 1

Main goals:

1. Design the data structure to represent the boost::real numbers.
2. Implement the general structure of the approximation intervals iterator (explained later).
3. Implement the explicit and operation number structures(explained later).
4. Implement the += and < operators for operation numbers (explained later)

Design the data structure to represent the boost::real::numbers

The first stage was very experimental, we first focus on researching about computable numbers and how to represent their intervals. It is worth to notice the fact that the bibliography regarding range arithmetic remains in an abstract mathematical dimension where we can refer to the error and describe the numbers intervals in terms of their approximation and the generated error. The problem with this in the practice is that the error and the aproximation are actually unknown, and because of this, we need to represent numbers intervals using another mechanism to work with the interval boundaries instead.

Note: To facilitate the research and the architecture design stage, we have used a decimal representation of the number and thus, all the digits are between 0 and 9. This decimal representation is sub-optimal in terms of performance and must modified to represent numbers in base 2³² so the full integer numbers are used.

In this stage, we have struggled for a while. We first tried to implement a structure that avoids the full boundaries recalculation each time the precision is increased. We decided to represent a number interval as the number truncation at a certain point for the interval lower boundary and the next available number using the truncated amount of digits for the interval upper boundary. This representation worked very well for those numbers explicitly described with no nines in their digits, as for example 1.0345. In this example, the worst approximation interval is [1, 2], the next approximation interval is [1.0, 1.1], the next approximation is [1.03, 1.04] and so on. For these cases, each time we iterate through the number intervals, we just need to:

1. Replace the last number of the upper boundary by the last number of the lowest boundary.
2. Append the next digit of the represented number to the lower boundary.
3. append the next digit of the represented number to the upper boundary and increase it by one.

But, what happens is the number is 1.0395? there is a moment when we have the interval approximation [1.03, 1.04] where if we follow the above algorithm we would generate carry when calculating the upper boundary and we will end up with [1.039, 1.03'10']. As we can see, this is not a valid number since 10 is not a symbol of the base 10. At this moment, we tried to get a method to avoid recalculating the entire upper boundary, but after several tries, we and up by figuring out that it is not possible. A very simple demonstration of this is as follows:

Consider a number a = 1.0…09…90…01 with an arbitrary amount of zeros at the beginning, after that an arbitrary amount of nines, then an arbitrary amount of zeros and finally a 1.

Then, we have the next three scenarios:

1. Any truncation made before the nines will not need to recalculate more than the previous digit and will return the interval [1.0…0, 1.0…1].
2. Any truncation after the first nine and before the last nine will return the interval [1.0…09…9, 1.0…10…0].
3. Any truncation made after the last nine and before the last 1 will return the interval [1.0…09…90…0, 1.0…09…90…1].

As we can see, in the second case, all the nines of the original numbers are replaced by zeros in the upper boundary while they remain as nines in the third case. Then, when passing from the second case to the third case, we are forced to recalculate all those digits. Because the original number could have an arbitrary number of consecutive nines, the boundary recalculation could imply recalculating an arbitrary number of digits. Thus, we cannot ensure a less than a linear complexity to recalculate the upper bound of a simple number.

After determining this result, we changed our algorithm to iterate through a number intervals backtracking in every iteration so the full upper bound is correctly recalculated.

Implement the explicit and operation number structures

In one hand a boost::real number can be constructed in many ways, for example, a number can be a literal, in which case, the number is written down in the code somewhere. Also a number can be created from a function (for example for the irrational numbers) or from the result of an operation between other two boost::real numbers. In the other hand, we want to represent a number as a set (finite or infinite) of intervals that approximate the number. Because we cannot calculate all the intervals (they can be infinite intervals) we construct the intervals each time the iterator is increased. If the number is a literal or an algorithm, then we construct the approximation intervals by simply truncating the numbers digits for the lower boundary and calculating the next representable number with that precision for the upper boundary as we have already explained, but if we have a number that comes from an operation of other two numbers, then we need to use the range arithmetic explained in [1] to calculate the interval of the operation. Because of this, each time an operation is calculated, we need a structure where to store copies of the operands so we can calculate the operation interval using its operands intervals. We need to keep copies of the operands because each time we want to iterate the operation, we first have to iterate the operands intervals and recalculate the new interval of the operation from the new intervals of the operands.

We have designed a binary tree structure, where the internal nodes represent the operations and the leaves are the final numbers that don’t come from an operation. We define then a boost::real number as one of the next three possibilities:

1. A real_explicit number, which is composed by:
   a. A std::vector<int> of digits.
   b. A bool specifying whether the number is positive (true) or negative (false).
   c. A decimal-point displacement (called exponent).
2. A real_algorithm number, which is composed by:
   a. A function or lambda function that returns the n-th digit of the number.
   b. A bool specifying whether the number is positive (true) or negative (false).
   c. A decimal-point displacement (called exponent).
3. A real_operation number which is an operation composed by
   a. Two boost::real pointers to its operands.
   b. An enumerate indicating the operation (+,- or *)

The Figure 5 shows an example of the binary tree representation of an operation.

Figure 5: A boost::real number resulting from a complex operations is represented by a binary tree.

The designed structure for the operations determines a recursive structure where several different instances can represent the same real number. For example if we consider the literal 2.0, then it is a real_explicit instance, but we can also consider 1.0 + 1.0 which is a real_operation with two literal children, or we can even consider (0.5 + 0.5) + (0.5 + 0.5) which is a real_operation with two real_operation as children.

Because so far we weren’t completely sure if the designed structure was good enough or if we would need to modify one more time the implemented algorithms, we decided to create a single C++ class with all the required attribute to represent any one of the three kinds of boost::real numbers (explicit, algorithm or operation) and refactor the code in the GSoC 2018 second stage once we got confident about our architecture.

Implement the general structure of the approximation intervals iterator

At this point of the project we had a structure of a binary tree that represents boost::real numbers, but as we have mentioned several times already, given a boost::real number we want to generate its approximation intervals from lower to higher precision so we can interact with these abstract numbers. For example, we would like to compare numbers, and because numbers can have an infinite amount of digits, we could not be able to compare them except by their intervals.

For this purpose, we worked on a const_precision_iterator that iterates the number set of approximation intervals. The iterator must start with the lower precision interval and increase the precision in each iteration. Once we have the iterator, numbers can interact using them. For example, to compare two numbers, we can first create iterators of those numbers and iterate their approximation intervals until they intervals stop overlapping and we are able to determine which number interval is lower than the other. As we will see later on, deciding if a real number is lower than another is not decidable, meaning that we can make programs that decide if one number is lower than the other if and only if both numbers are not equal.

Note: The idea is not that the boost::real user manipulates the iterator. Instead, these iterators will be internally used by the operations. For example, when trying to determine if a number x is lower than another number y, the developer just write “x < y”, and inside the < operator, the numbers iterators will be used to determine the value of the operator.

Now, we will introduce the boost::real::const_precision_iterator data type that we have designed and implemented in the first stage of the GSoC 18. The iterator is a const iterator because it does not modify the number that is iterating, thus, we are able to implement the < operator which is a const operator in a safe mode.

The next code show the const_precision_iterator attributes:

// Internal number to iterate
real const* _real_ptr = nullptr;
// Explicit or algorithmic number iterator
boost::real::interval interval;
// If the number is a composition, the const_precision_iterator uses the operand iterators
boost::real::const_precision_iterator* _lhs_it_ptr = nullptr;
boost::real::const_precision_iterator* _rhs_it_ptr = nullptr;

A pointer to const is used to ensure the iterator will not change the number that is pointing. Also a boost::real::interval structure is used to store the current approximation interval of the iterator. An interval is a struct with a lower and upper bound. Finally, if the pointed real number is an operation, then we have the iterators of its operands, so we can use those iterators to recursively calculate its operand intervals and from those intervals, it calculates the current number interval in different ways depending on the performed operation (+- or *).

A boost::real::interval is defined as follows:

struct interval {
    boost::real::boundary lower_bound;
    boost::real::boundary upper_bound;
}

And a Boost::real::boundary is defined as follows:

struct boundary {
    std::vector<int> digits = {};
    int exponent = 0;
    bool positive = true;
}

A boundary has indeed, the same structure than a boost::real number when it is an explicit number.

Then, the algorithm to iterate over the const_precision_iterator of a boost::real number has a simple general schema shown in the next pseudo-code:

operator++():
    if is a real_explicit or algorithm number:
        1. Add one more digit to the interval lower bound truncation
        2. Recalculate the interval upper bound
    
    else
        1. increase operands iterators.
        2. calculate new interval boundaries depending the operation

Implement the += and < operators for operation numbers

Now that we got the general structure of the boost::real numbers representation and its approximation interval iterator, we started constructing the mathematical operators which are the final goals of working with real numbers.

Because a boost::real number can be one of real_explicit, real_algorithm or real_operation, when we apply a not const operator to a number, the number can be converted from one type to another. For example, if we have the real_explicit 3.0 and we use the += operator with another number:

1. boost::real x = "3.0";
2. boost::real y = "2.0";
3. x += y

We have that after the line number 1, x is a real_explicit, but after the line number 3, x is converted, because now it is a real_operation where y and itself are the operands. To achieve this, and because the new real_operation number must have valid pointers to its operands without depending on whether they are externally destroyed or not, it must first make deep-copies of the operands (the right side operand and itself) and then convert the number into an operation correctly setting the _kind and _operation attribute to declare which operation it represents.

Figure 6: Graph representation of the operation +=. A real_explicit number is transformed to a real_operation binary tree with two children as the leaves (copies of the operands) and the operation kind as KIND::OPERATION. Also the operation is set to OPERATION::ADDITION.

The code is as follows:

real& operator+=(const real& other) {
    this->_lhs_ptr = new real(*this);
    this->_rhs_ptr = new real(other);
    this->_kind = KIND::OPERATION;
    this->_operation = OPERATION::ADDITION;
    return *this;
}
real& operator-=(const real& other) {
    this->_lhs_ptr = new real(*this);
    this->_rhs_ptr = new real(other);
    this->_kind = KIND::OPERATION;
    this->_operation = OPERATION::SUBTRACTION;
    return *this;
}
real& operator*=(const real& other) {
    this->_lhs_ptr = new real(*this);
    this->_rhs_ptr = new real(other);
    this->_kind = KIND::OPERATION;
    this->_operation = OPERATION::MULTIPLICATION;
    return *this;
}

Until here, we have generated the number representations, these representations are abstract because they do not explicitly represent the numbers but instead it represents the whole operation binary tree from where the number results. Then, in order to use these numbers, we need to implement the approximation intervals algorithm, which is different for each operation.

The interval calculation for the += operator
Once, the boost::real number is converted into an addition operator, we calculate the approximation interval returned by the precision iterator algorithm from the boundaries of its operands as follows:

void calculate_addition_interval_boundaries() {
    
    boost::real::helper::add_boundaries(
            this->_lhs_it_ptr->interval.lower_bound,
            this->_rhs_it_ptr->interval.lower_bound,
            this->interval.lower_bound
    );
    boost::real::helper::add_boundaries(
            this->_lhs_it_ptr->interval.upper_bound,
            this->_rhs_it_ptr->interval.upper_bound,
            this->interval.upper_bound
    );
}

where boost::real::helper::add_boundaries is a function that implements the addition between two boost::real::boundary variables, considering the exponent to align both boundaries digits, the full code of this helper method does not clarify the idea of the algorithm and can be seen here.

The logic behind this algorithm is simple: if we have the interval of two numbers x1 ≤ x ≤ x2 and y1 ≤ y ≤ y2, then, because the addition is an increasing monotone function, we now that x1 + y1 ≤ x + y ≤ x2 + y2. Then, the interval of the number x + y can be set as [x1 + y1, x2 + y2] without knowing what is the exact value of x and y. This is why, in the code, we calculate the addition of the lower boundaries as the interval lower boundary and the addition of the upper boundaries as the interval upper boundary.

The < operator
The < operation attempt to determine if the left side operator is less than the right side operator. Because we are dealing with intervals, we need to:

1. Determine what does it mean for an interval to be lower than another interval in a way that is consistent with the definition for the real numbers represented by the intervals.
2. Determine which precision are we going to consider.

In interval arithmetic, an interval x is lower than another interval y if and only if the upper boundary of x is lower than the lower boundary of y. In other words, x < y if and only if, for every number in x it is lower than any number in y.

If a number x is not equal to the number y, because iterators increases the precision by making each time the intervals smaller, eventually the intervals will stop overlapping. Then, we need to iterate until we reach a point where intervals stop overlapping.

As it is better explained in [1], determining if a real number is lower than another real number is not decidable, meaning that if the compared numbers are not equal, we can decide, but if the numbers are equals, we can’t decide. The major problem with this is that we can never know when to stop trying, and if we stop trying we would never know whether we stop trying to soon or if both numbers are equals. Because of this, we need to determine a maximum precision to stop if the number intervals still overlap at that precision. When the maximum allowed precision is reached, the operation will fail and throws a precision exception. When a precision exception occurs, the user can increase the maximum allowed precision or abort the operation. Also, the maximum precision can be modified in run-time by a custom logic defined to handle the maximum precision.

Figure 7: precision iteration example, x is lower than y, but for the first two iteration over the approximation intervals, the intervals of x and y still overlap and the decision can’t be made. In the third iteration, the number intervals are small enough to stop overlapping and we are able to compare x with y.

The operator code is as follows:

bool operator<(const real& other) const {
    auto this_it = this->cbegin();
    auto other_it = other.cbegin();

    int current_precision = std::max(
        this->max_precision(),
        other.max_precision()
    );
    for (int p = 0; p < current_precision; ++p) {
        // Get more precision
        ++this_it;
        ++other_it;

        if (this_it.interval < other_it.interval) {
            return true;
        }

        if (other_it.interval < this_it.interval) {
            return false;
        }
    }

    // If the precision is reached and the number intervals 
    // still overlap, then we cannot determine if they are equals 
    // or other is less than this using this maximum precision, 
    // so we throw an error.
    throw boost::real::precision_exception();
}

GSoC stage 2

Main goals:

1. Implement the -= operator interval calculation.
2. Implement the *= operator interval calculation.
3. Implement the +, - , * and = operators.
4. Implement the real_algorithm number type.
5. Code refactor to separate the boost::real class in three (boost::real_explicit, boost::real_algorithmic and boost::real).

Implement the -= operator interval calculation

For the interval of a subtraction between two intervals consider x1 ≤ x ≤ x2 and y1 ≤ y ≤ y2. Then, because the subtraction is a decreasing monotone function, we know that x1 - y2 ≤ x - y ≤ x2 - y1:

void calculate_subtraction_interval_boundaries() {
    
    boost::real::helper::subtract_boundaries(
            this->_lhs_it_ptr->interval.lower_bound,
            this->_rhs_it_ptr->interval.upper_bound,
            this->interval.lower_bound
    );
    boost::real::helper::subtract_boundaries(
            this->_lhs_it_ptr->interval.upper_bound,
            this->_rhs_it_ptr->interval.lower_bound,
            this->interval.upper_bound
    );
}

Implement the *= operator interval calculation

Once again, to determine the boundaries of the interval of a multiplication we need as for the other operations to calculate an interval containing every possible value resulting from multiplying numbers of the operands intervals. i.e. if we have x1 ≤ x ≤ x2 and y1 ≤ y ≤ y2, then we want to calculate z1 and z2 from x1, x2, y1 and y2 to satisfy z1 ≤ x*y ≤ z2. The problem is that multiplication is not a monotone function and we must consider the sign multiplication:

• (+) * (+) = (+)
• (+) * (-) = (-)
• (-) * (+) = (-)
• (-) * (-) = (+)

This is commonly easy to remember (equal sign results in positive numbers and different signs results in negative numbers) but; how do we calculate the boundaries of the multiplication number from the boundaries of its operands?

We have separate the interval calculation in the next cases:

1. Both intervals are fully contained in the positive number line.
2. Both intervals are fully contained in the negative number line.
3. One interval is fully contained in the negative number line and the other interval is fully contained in the positive number line.
4. There are boundaries with different sign.

Positive intervals:
If both intervals are positive we know that x1 ≤ x ≤ x2 and y1 ≤ y ≤ y2 implies that x1*y1 ≤ x*y ≤ x2*y2 because in the positive number line, the multiplication is monotone.

Negative intervals:
If both intervals are negative, then we know that x1 ≤ x ≤ x2 and y1 ≤ y ≤ y2 implies x2*y2 ≤ x*y ≤ x1*y1. This is because multiplication of negative numbers results in a positive number, then:

• x1 ≤ 0 and y1 ≤ 0 implies that x1*y1 = |x1|*|y1|
• x2 ≤ 0 and y2 ≤ 0 implies that x2*y2 = |x2|*|y2|
• x1 ≤ 0, x2 ≤ 0 and x1 ≤ x2 implies |x1| ≥ |x2|
• y1 ≤ 0, y2 ≤ 0 and y1 ≤ y2 implies |y1| ≥ |y2|

As a result of the above implication we have that:

• |x1|*|y1| ≥ |x1|*|y2|
• |x1|*|y1| ≥ |x2|*|y1|
• |x1|*|y1| ≥ |x2|*|y2|
• |x2|*|y2| ≤ |x1|*|y2|
• |x2|*|y2| ≤ |x2|*|y1|
• |x2|*|y2| ≤ |x1|*|y1|

Positive vs Negative interval:
For the case of a negative interval [x1, x2] multiplied by a positive interval [y1, y2] (An analogous reasoning works for the opposite case) we have that:

• x1 ≤ x ≤ x2 and y1 ≤ y ≤ y2 implies x*y ≤ 0.
• x1 ≤ 0, x2 ≤ 0 and x1 ≤ x2 implies |x1| ≥ |x2|
• y1 ≥ 0, y2 ≥ 0 and y1 ≤ y2 implies |y1| ≤|y2|

As a result of the above implications we have that:

• |x1|*|y2| ≥|x1|*|y1|
• |x1|*|y2| ≥ |x2|*|y1|
• |x1|*|y2| ≥ |x2|*|y2|
• |x2|*|y1| ≤ |x1|*|y1|
• |x2|*|y1| ≤ |x2|*|y1|
• |x2|*|y1| ≤ |x2|*|y2|
• x2*y1 ≤ 0
• x1*y2 ≤ 0

And the thus, x1 ≤ x ≤ x2 and y1 ≤ y ≤ y2 imples that a x2*y2 ≤ x*y ≤ x1*y2.

Boundaries with different sign:
Finally, if we have two boundaries with a different sign, we need to calculate the four possibilities x1*y1,x1*y2, x2*y1 and x2*y2 and test for the minimum and maximum of those numbers, we know that a*b must be between those boundaries.

The interval calculation code for the multiplication is as follows:

void calculate_multiplication_interval_boundaries() {
    
    // positive() return true if and only if the full 
    // interval is positive
    bool lhs_positive = this->_lhs_it_ptr->interval.positive();
    bool rhs_positive = this->_rhs_it_ptr->interval.positive();
    // negative() return true if and only if the full interval 
    // is negative
    bool lhs_negative = this->_lhs_it_ptr->interval.negative();
    bool rhs_negative = this->_rhs_it_ptr->interval.negative();
   if (lhs_positive && rhs_positive) { 
        
        // positive interval vs positive interval 
        boost::real::helper::multiply_boundaries(
                this->_lhs_it_ptr->interval.lower_bound,
                this->_rhs_it_ptr->interval.lower_bound,
                this->interval.lower_bound
        );
        boost::real::helper::multiply_boundaries(
                this->_lhs_it_ptr->interval.upper_bound,
                this->_rhs_it_ptr->interval.upper_bound,
                this->interval.upper_bound
        );
    } else if (lhs_negative && rhs_negative) { 
        // negative interval vs negative interval
        boost::real::helper::multiply_boundaries(
                this->_lhs_it_ptr->interval.upper_bound,
                this->_rhs_it_ptr->interval.upper_bound,
                this->interval.lower_bound
        );
        boost::real::helper::multiply_boundaries(
                this->_lhs_it_ptr->interval.lower_bound,
                this->_rhs_it_ptr->interval.lower_bound,
                this->interval.upper_bound
        );
    } else if (lhs_negative && rhs_positive) {
     
        // negative interval vs positive interval
        boost::real::helper::multiply_boundaries(
                this->_lhs_it_ptr->interval.lower_bound,
                this->_rhs_it_ptr->interval.upper_bound,
                this->interval.lower_bound
        );
        boost::real::helper::multiply_boundaries(
                this->_lhs_it_ptr->interval.upper_bound,
                this->_rhs_it_ptr->interval.lower_bound,
                this->interval.upper_bound
        );
    
    } else if (lhs_positive && rhs_negative) {
   
        // positive interval vs negative interval
        boost::real::helper::multiply_boundaries(
                this->_lhs_it_ptr->interval.upper_bound,
                this->_rhs_it_ptr->interval.lower_bound,
                this->interval.lower_bound
        );
        boost::real::helper::multiply_boundaries(
                this->_lhs_it_ptr->interval.lower_bound,
                this->_rhs_it_ptr->interval.upper_bound,
                this->interval.upper_bound
        );
    } else {
        // one interval is around zero and then all 
        // possible combinations are tested
        boundary current_boundary;
        // Lower * Lower
        boost::real::helper::multiply_boundaries(
                this->_lhs_it_ptr->interval.lower_bound,
                this->_rhs_it_ptr->interval.lower_bound,
                current_boundary
        );
        this->interval.lower_bound = current_boundary;
        this->interval.upper_bound = current_boundary;
        // Upper * upper
        boost::real::helper::multiply_boundaries(
                this->_lhs_it_ptr->interval.upper_bound,
                this->_rhs_it_ptr->interval.upper_bound,
                current_boundary
        );
        if (current_boundary < this->interval.lower_bound) {
            this->interval.lower_bound = current_boundary;
        }
        if (this->interval.upper_bound < current_boundary) {
            this->interval.upper_bound = current_boundary;
        }
        // Lower * upper
        boost::real::helper::multiply_boundaries(
                this->_lhs_it_ptr->interval.lower_bound,
                this->_rhs_it_ptr->interval.upper_bound,
                current_boundary
        );
        if (current_boundary < this->interval.lower_bound) {
            this->interval.lower_bound = current_boundary;
        }
        if (this->interval.upper_bound < current_boundary) {
            this->interval.upper_bound = current_boundary;
        }
        // Upper * lower
        boost::real::helper::multiply_boundaries(
                this->_lhs_it_ptr->interval.upper_bound,
                this->_rhs_it_ptr->interval.lower_bound,
                current_boundary
        );
        if (current_boundary < this->interval.lower_bound) {
            this->interval.lower_bound = current_boundary;
        }
        if (this->interval.upper_bound < current_boundary) {
            this->interval.upper_bound = current_boundary;
        }
    }
}

Implement the +, - , * and = operators

Once we have the operators +=, -= and *=, we can easily implement the const operators +, — and *. The main difference between these operators is that they don’t modify the right side operand. Instead, they create a new boost::real number containing the result of applying the operands:

real operator+(const real& other) const {
    real result = *this;
    result += other;
    return result;
}
real operator-(const real& other) const {
    real result = *this;
    result -= other;
    return result;
}
real operator*(const real& other) const {
    real result = *this;
    result *= other;
    return result;
}

For the = operator, we need to consider the fact that if the boost::real number to copy is an operation, then it contains pointers to its operands. If we want to avoid aliasing we need to implement deep-copies by first making new copies of the operands and store pointers to those new copies.

real& operator=(const real& other) {
    this->_kind = other._kind;
    this->_explicit_number = other._explicit_number;
    this->_operation = other._operation;
    
    if (other._lhs_ptr != nullptr) {
        this->_lhs_ptr = new real(*other._lhs_ptr);
    }
    if (other._rhs_ptr != nullptr) {
        this->_rhs_ptr = new real(*other._rhs_ptr);
    }
    return *this;
}

Copying a boost::real number operation requires to recursively copy its operands. Because of this, the copy operator algorithm has a linear complexity over the number of operations from which the number was constructed. Therefore, copying more complex numbers is less efficient than copying simpler numbers even if the numbers represent the same value.

Code refactor to separate the boost::real class in three (boost::real_explicit, boost::real_algorithmic and boost::real)

As we pointed out before, a boost::real number is a binary tree that represents the operation from where the number results and a number have infinite representations. In order to obtain a more clear code, we decided to separate the class into three different C++ classes we have reached the last and higher possible precision intervalto represent each different kind of boost::real number (real_explicit, real_algorithm and real_operation).

Boost::real_explicit:
This class stores a not recursive explicit number, meaning that the number digits are explicitly defined in a vector. The class attributes are the same as the boost::real::boundary.

// Number representation as a vector of digits with an 
// integer part (_exponent) and a sign +/- (_positive).
std::vector<int> _digits;
int _exponent;
bool _positive;

The boost::real_explicit::const_precision_iterator calculates the lower boundary as a truncation of the original number and the upper boundary as the next valid number using that number of digits as explained before. The code is as follows:

void operator++() {

    // If the explicit number full precision has 
    // been already reached,
    // then just add zeros at the boundaries ends
    if (this->_n >= (int)this->_real_ptr->_digits.size()) {
        this->_n++;
        return;
    }

    // If the number is negative, boundaries are 
    // interpreted as mirrored: First, the operation 
    // is made as positive, and after boundary calculation
    // bounds are swapped to come back to the negative 
    // representation.
    this->check_and_swap_boundaries();

    // If the explicit number just reaches the full precision
    // then set both boundaries equals
    if (this->_n == (int)this->_real_ptr->_digits.size() - 1) {

        this->interval
            .lower_bound
            .push_back(this->_real_ptr->_digits[this->_n]);
        this->interval
            .upper_bound = this->range.lower_bound;


    } else {

        // If the explicit number didn't reaches the full precision
        // then the number interval is defined by truncation.
        this->range
            .lower_bound
            .push_back(this->_real_ptr->_digits[this->_n]);

        this->range.upper_bound.clear();
        this->range
            .upper_bound
            .digits
            .resize(this->range.lower_bound.size());
        // Calculate the upper bound        
        int carry = 1;
        for (int i = (int)this->interval.lower_bound.size() - 1; 
             i >= 0; 
             --i) {
            
            if (this->interval.lower_bound[i] + carry == 10) {
                this->interval.upper_bound[i] = 0;
            } else {
                int digit = this->range.lower_bound[i] + carry;
                this->interval.upper_bound[i] = digit;                
                carry = 0;
            }
        }

        // If we still have a carry in the most left digit, 
        // then we add a new digit to the number.
        if (carry > 0) {
            int new_digit = this->range.lower_bound.exponent + 1;            
            this->interval.upper_bound.push_front(carry);
            this->interval.upper_bound.exponent = new_digit;
        
        } else {
            int exponent = this->range.lower_bound.exponent;
            this->interval.upper_bound.exponent = exponent;
        }
    }

    // If the number is negative, then the boundaries are mirrored
    // and must comeback to their original order.    
    this->check_and_swap_boundaries();
    this->_n++;
}

The upper boundary recalculation complexity is linear until the precision is higher than the number digit amount. After the precision is higher than the number of digits, the interval lower and upper boundaries remain equals and we have reached the perfectly accurate representation of the number.

Boost::real_algorithm:
This class is exactly as the real_explicit number with the difference that instead of having a vector with the digit explicitly defined, we have a pointer to a function or a lambda function (_get_nth_digit) that for each unsigned int “n” the function returns the value of the n-th digit.

// Instead of using:
this->_real_ptr->_digits[this->_n]
// We use:
this->_real_ptr->_get_nth_digit(this->_n)

Boost::real (real_operation):
This is the main class that represents all the boost::real number operation binary tree. At this point, we have tried to separate the real_operation from the main class as we did with the boost::real_explicit and boost::real_algorithmic, but there is a problem; A real operation is actually an internal or root node in the binary tree and it has two boost::real numbers pointers to its operands, those numbers can be of any type (real_explicit, real_algorithmic or real_operation), then this class has a recursive definition and we need the boost::real definition to define the class. At the same time, because the boost::real is the main class, then it must be able to store any kind of numbers and therefore, it needs the boost::real_operation definition as well. Because of this cyclic dependency, we were not able to separate the real operation class from the main boost::real and after some trials, we ended up by defining the real_operation in the main boost::real class.

Even if we prefer to have a separated boost::real_operation class to enforce the semantic, having the nodes of the tree as the main boost::real class is actually the correct pattern when dealing with graph structures as trees and lists.

A boost::real class is composed by:

// Available operations
enum class OPERATION {ADDITION, SUBTRACT, MULTIPLICATION, NONE};
enum class KIND {EXPLICIT, OPERATION, ALGORITHM};

KIND _kind;

// Explicit number
real_explicit _explicit_number;

// Algorithmic number
real_algorithm _algorithmic_number;

// Composed number
OPERATION _operation = OPERATION::NONE;
real* _lhs_ptr = nullptr;
real* _rhs_ptr = nullptr;

The two enumerate classes allows us to define which type of number the boost::real is representing.

• KIND specifies if the number is a real_explicit, a real_algorithm or a real_operation number.
• OPERATION specifies which operation (+, - or *) does the number represents.

In case the number is an operation, the class has two recursive pointers to boost::real numbers and it specifies the operation (ADDITION, SUBTRACTION or MULTIPLICATION).

If the represented number is an explicit or algorithmic number, then all the methods will forward the operation to the corresponding method of the respective classes. If the number is an operation, then, the class implements the recursive algorithms we have described for each operation.

GSoC stage 3

Main goals:

1. Explore the /= operator.
2. Code refactor to correctly use the exponent as a logical alignment of the number digits (explained later).
3. Implement more constructors, overload the << operator and add a method to set the maximum precision.
4. Implement operators == and >.
5. Implement an irrational number as example that comes with the library.
6. Implement tests.
7. Document the new library.

Explore the /= operator

We have left the /= operator as an experimental stage because of the next issues:

1. Dividing two vectors is more complex than for the other operations.
2. The division of two finite numbers can result in an infinite number (consider 1/3 = 0.33…)
3. The calculation of the exponent is not simple unless we consider the entire integer part of the result which could be larger than the maximum precision.

We were able to address the first two items: vectors division algorithm was implemented specifying the maximum amount of considered digits of the result and performing a truncation if that maximum is reached. This approach has a small issue; If we are calculating an upper bound, performing a truncation will give us a lower number than the real boundary and the represented number could not be represented. To solve this issue, we must calculate the next available number with that number of digits as we did for the real_explicit number.

The problem we couldn’t solve yet for the division algorithm is how to calculate the exponent of the division algorithm if we still are truncating the number before reaching the entire integer part. Consider the division of 1000 / 2, the result is 500 which is 0.5x1⁰³, the exponent is three, but what happens if we are just considering 2 digits of the result? Then if we follow the well-known division algorithm we have that:

1. 10 > 2, then 2 * 5 = 10 and we set the result first digit as 5.
2. we have the remaining of 10 - 2*5 = 0, and because 0 < 2 we add a zero in the result and we got our two digits 50, but at this point, we need to know if the exponent is 2 or another higher number.

We have then, explore the possibilities of the division algorithm and we are close to a solution, but this will be implemented after the GSoC 18 and will be committed to the library separately from the rest of the project so we can implement a solid division algorithm.

The vector division algorithm that we have already implemented is stored here.

Code refactor to correctly use the exponent as a logical alignment of the number digits

Until this point in order to simplify the algorithms each time two vector where added, subtracted or multiplied, we first created two new vectors and align those new vectors so the vector indexes represent the same units; Let’s see the next example:

• 10.25 has the vector of digits <1,0,2,5> and the exponent 2.
• 0.041 has the vector <4,1> and the exponent -1.

If we add these two numbers, we get the result 10.25 + 0.041 = 10.291:

10.250
00.041
— — —
10.291

as we can see, the last element of the 10.25 vector (the number “5”) is aligned with the first element of the 0.041 vector (the number “4”). Because of this, we first transformed the vectors to have the same length and represent the same units as follows:

• <1,0,2,5> is transformed to <1,0,2,5,0>
• <4,1> is transformed to <0,0,0,4,1>

Now, the 5 and the 4 are in the same index (they are both the four elements o the arrays). This solution was good to move forward when we just started designing the new library from scratch and we needed to design the basic structure of the data type. But it is very inefficient in both, time and memory complexity.

The solution is to consider a universal logical position where the first fractional number is at position zero and positions increases in the direction to lower digits and decreases in direction of higher digits, then, all the digits of the integer part of the number will have a negative logical position. To get the real position of the digit in the vector, we need to add the exponent of the number to the logical position, and if the real position is not a valid index the digit is a zero.

In the case of the number 0.041, which has the vector <4,1> and exponent -1 we can see that the logical positions -1 and 0 goes to -2 and -1 respectively, and because they are not valid position in the vector <1,4> then we know that those digits are both zeros.

Implement more constructors, overload the << operator and add a method to set the maximum precision

Considering the fact that in floating point representations, not all the numbers are valid, each time we initialize a double variable with a literal number we could face an approximation error, then, it is not safe to initialize a boost::real from a double number and we need a constructors that remain intuitive and without losing precision, as for example, constructing from a string:

boost::real a = “-1.3404”

The example should return a boost::real_explicit number with the _digits vector <1,3,4,0,4>, the _exponent as 1 and the _positive as false.

For this purpose, we have created a constructor that takes a string representing a number in base 10 and generate a real_explicit number representation from it. We have also overloaded the = operator for std::strings.

We also implemented different constructors that use initialization lists for those cases where performance is important to avoid the conversion from string to std::vector<int>. All the available constructors are explained in the GitHub project README.md file.

Finally, we have set a static member to the boost::real::real and boost::real::real_algorithm classes to set the maximum precision. Because of this, the developer must initialize the maximum precision with an unsigned int before being able to use the library as shown in the following example:

unsigned int boost::real::real::maximum_precision = 10;
unsigned int boost::real::real_algorithm::maximum_precision = 10;

If in any moment of the program execution the maximum precision need to change, it only requires to change the value of these static members.

Because we could need to use different maximum precision for different instances, we can use the instance method set_maximum_precision(unsigned int) to set a particular maximum precision for the instance and the static members maximum precision are used as the default precision when the instance has not set a particular maximum precision. Note: Setting a maximum precision of zero causes the instance to use the default static maximum precision.

Implement operators == and >

The > operator was implemented analogous to the < operator, the difference is that an interval is greater than another interval if and only if the lower boundary of the interval is greater than the upper boundary of the other interval.

For the == operator we must either reach the full number precision or reach a precision where the intervals stop overlapping. Then, we iterate both operands until one of the next events happen:

1. The operator got the full operand number precisions and we can exactly compare both numbers.
2. The operator reaches a precision where the interval stop overlapping and determine that they are not equals.
3. The operator reaches the maximum precision and has not yet reached the full precision of both operands. In this case, the operator throws a precision exception to indicate that for the set maximum precision it cannot determine if the numbers are equals or not

Implement an irrational number as example that comes with the library

We have decided to implement a simple irrational number to come with the library as an example for future users. We also look forward to adding some important irrational numbers as π to come with the library because they are commonly used.

For this purpose, we have created the boost::real::irrational namespace and placed the Champernowne irrational number for binary representation. The Champernowne number is defined by successively concatenating the integers numbers in a determined base. For example, in the case of the binary base we have the integers numbers:

0 1 10 11 100 101 110 111 1000 1001 1010 1011 1100 1101 1110 1111

Then, the Champernowne number for binary representation is as following:

0.1101110010111011110001001101010111100110111101111

The function that returns the n-th digit of the Champernowne number for binary representation starts with the binary number 0 increase the number while an iterator iterate through the number digits until it gets to the n-th digit and it returns that digit.

Implement tests

Testing the iterators:
Testing is an entire world itself, ideally, we would make white and black box testing, and for the black box testing, we would need a separated person not familiarized with the code. Because this is not the case for a Google Summer of Code we have implemented tests in the middle between white and black box testing. We used some information of how does the data type works to make tests of the border cases and other general cases and cover the most important part of the algorithm.

We decided to use the Catch2 library that provides a single file header-only include because it provides us with enough power to test the library and it also removes dependencies for the future users of the library. Also, we decided to use Travis-CI for continuous integration, and Catch2 worked very well with Travis-CI.

We have separated tests in unit tests for the boost::real_explicit and boost::real_algorithm data types and integration test for the boost::real.

For the unit test, we made tests for the iterator using combinations of:

• The number sign: + and -
• The number exponent: -3, -2, -1, 0, 1, 2 and 3.
• With carrying and without carrying: We used number witch one or multiple 9s as digits to ensure the operation will generate carry and numbers with all digits lower than 5 to never have carry.

For the integration test, we designed two series of tests for the iterators:

1. boost::real numbers of a single operation for different combinations of numbers.
2. boost::real numbers that combines operations.

Test the base cases of the recursive structure and the first level gives us confidence in the basis of all the recursive structure. The structure next levels are abstracted from their position within the binary tree and have the same logic as the second label.

For the tests of a single operation, we considered combinations of:

• Operand sign: (+,+) (+,-) (-,+) and (-,-)
• Operator: +, - and *.
• Operand type: real_algorithm and real_explicit.
• With carrying and without carrying.

For the tests of combined operation we have considered the next properties:

• Operator combinations: (+,+) (+,-) (+,*) (-,+) (-,-) (-,*) (*,+) (*,-)(*,*)
• With carrying and without carrying.
• (x op1 y) op1 z vs x op1 (y op2 z), where op1 and op2 are operators +, - and *.

We would rather prefer to test the combinations with and without carrying than testing different operand sign or testing the operand types (real_explicit or real_algorithm). The reason for this is that the number sign or type don’t expand throw the recursive structure and the algorithm doesn’t significantly change for those case, in the other hand, the carry is a main part of the interval calculation algorithms and carries error an easily be propagated through the operations.

Test the operators:
Because the final goal of a numeric number is to test the operators that allows us to work we the number, we have defined a battery of tests for the operators <, == and > considering the next attributes:

• The numbers sort (x = y, x < y, x > y).
• The numbers intervals still overlap after the maximum precision is reached vs the numbers intervals stop overlapping before the maximum precision is reached.
• The number representation (Explicit, Algorithm, Addition (x = x1 + x2), subtraction (x = x1 - x2) and multiplication (x = x1*x2).

Document the new library

We have documented this first version of the library in three different places:

• README.md: The library introduction and reference manual with the API are documented in the README.md file.
• Doxygen: We made a documentation of the data structures and method specifications in Doxygen format within the code and exported as HTML. This documentation is more technical and it is intended for those developers interested in working in the library.
• This medium story: In this GSoC 18 project, the algorithms are important but the ideas behind those algorithms are more important. The entire project was a constant exploration of possibilities and abstractions. In this story, we hope to show the main ideas of the library and an overview of the creation process of a new library.

Future work

Apart from finishing the division algorithm we explored in this project we have some performance improvements we would like to do.

Change the boost::real number base from 10 to 2³² so we don’t waste all the number digit memory:
In the current version of the library we use integer numbers to represent the digits of decimal numbers, because of this, we are wasting lot of memory by only using 4 bits of the 32 bits of a integer number of the C++ int_32. To solve the memory problem of this issue, we must first implement a conversion algorithm to consume and print numbers in decimal representation (for usability purposes)but to internally represent numbers in base 2³². Also, we must update the addition subtraction and multiplication algorithm to work with numbers in base 2³².

Compress the binary tree when possible:
If we are working with a real_operation number composed by two operands and we are iterating over the number precision intervals and we reach the point where the lower boundary of the number is equal than the upper boundary of the number, then we know that the number is exactly equal to its boundaries, and we are able to replace the real_operation number by the new real_explicit representation. In order to implement this, we need a not const iterator, because iterations will have a side effect on the iterated number representation.

Use std::variant in the boost::real class:
The boost::real class is actually a class that either is a real_operation and has pointers to its two operands, or is one of real_explicit or real_algorithm and has pointers to boost::real_explicit and boost::real_algorithm variables respectively. This uses more memory than needed because it always has unused attributes. For example, when we have a real_explicit number, all the boost::real members that are not related with a real_explicit number are not used.

In order to improve the memory consumption, we want to move the structure to an std::variant schema that will correctly handle the memory to use the correct representation depending on the number kind (boost::real_explicit, boost::real_algorithm and boost::real).

Share number pointers, converting the binary tree into a DAG to reuse representations:
As explained in [2], if we have the next code

boost::real x = "1";
for (i=0; i<100; ++i) {
    x = x+x;
}

Then x is a number represented in a binary tree composed of 2¹⁰¹ - 1 nodes that are all equals. This performance problem is not easy to solve but can be addressed by using shared pointers and determine when a number is mutable or not in order to share the pointer with other numbers that want to use it. This improvement would convert the binary tree representation by a DAG (directed acyclic graph).

Implement Benchmarks:
The designed structure and algorithm have a higher temporal and memory complexity than using native number representation as the floating point and even if in the theory we know the complexity of the algorithms, we should implement benchmarks to check the complexity in the practice and be able to correctly profile and improve the performance in future releases.

References

1. Computable calculus / Oliver Aberth, San Diego : Academic Press, c2001
2. Lambov, B. (2007). RealLib: An efficient implementation of exact real arithmetic. Mathematical Structures in Computer Science, 17(1), 81-98.
3. Aberth, O., & Schaefer, M. J. (1992). Precise computation using range arithmetic, via C++. ACM Transactions on Mathematical Software (TOMS), 18(4), 481-491.
4. https://en.cppreference.com/w/cpp/concept/ForwardIterator
5. Damián Alberto Vicino. Improved time representation in discrete-event simulation. Other [cs.OH]. Université Nice Sophia Antipolis, 2015. English. .