Bringing Boost.Real to Review Ready State
Google Summer of Code 2020
Documentation:
Information:
Source Code: https://github.com/BoostGSoC20/Real
Introductory blog: https://medium.com/@laobelloli/boost-real-9e2dfbfbed5b(Created by Lao Belloli)
Mentors: Damian Vicino, Laouen Belloli
Partner in GSoC’20: Kishan Shukla
Project Completion: During GSoC’ 20, I completed all work that I got from my proposal to bring Boost.Real to review ready state. I implemented a new and optimized system for integers and rational numbers(integer_number class and real_rational class). I created Taylor expansions for various functions(exp, logarithm, sin, cos, tan, etc). Then used log and exp to make power function. I added several constants required for basic scientific calculations. After combining mine and Kishan Shukla’s work, I hope this library will be ready for review for adding it to boost libraries.
Commit History: https://github.com/BoostGSoC20/Real/commits?author=vikram2000b
Pull Requests: https://github.com/BoostGSoC20/Real/pulls?q=is%3Apr+author%3Avikram2000b
Introduction:
This library aims at providing a data-type that can represent all types of computable real numbers and provides the flexibility of performing arbitrary precision arithmetic, the user can always decide and control the amount precision to be used for a calculation thus controlling the computational power used by the machine.
A number is represented by a program/class x(k) that returns the intervals containing x. Those intervals are sorted from lower to higher precision, the precision of those intervals will increase as we keep on iterating for more precise results. And, we can always decide till the iterations that should be used for the calculation, eventually deciding the maximum precision used by the machine.
All calculations in this library are done using range arithmetic as every number is represented by an interval.
For example, x = 1.2356
In decimal base, the intervals will be something like this:
First iteration, x = [1, 2]
Second iteration, x = [1.2, 1.3]
Third iteration, x = [1.23, 1.24]
Fourth iteration, x = [1.235, 1.236]
Fifth and last iteration, x = [1.2356, 1.2357].
Maximum precision is when we can not produce more precise intervals which can represent our number. Here we can not iterate after the fifth iteration, as we have reached maximum precision. The maximum number of iterations can always be customized according to the power of the machine.
The foundational work was done in GSoC’18 by Lao Belloli and Damian Vicino, then some improvements were done in the library in GSoC’19 by Sagnik Dey and Kimberly Swanson. The link to the blog of everyone’s work is:
- Laouen Belloli (GSoC’18): https://medium.com/@laobelloli/boost-real-9e2dfbfbed5b
- Sagnik Dey (GSoC’19): https://sagnikdey92.github.io/GSoC
- Kimberly Swanson (GSoC’19): https://universenox.github.io/gsoc19_Final_Eval
In GSoC’20, I worked with Kishan Shukla. Link of Kishan Shukla’ work:
- Kishan Shukla (GSoC’20): https://kishanshukla-2307.github.io/posts/2020/08/gsoc-final-evaluation/
This library represents every number using intervals (according to precision) and then do all calculations accordingly. Initially, three data-types were implemented to represent a real number:
- Real Explicit: This is used when the number is a literal, number is written down in code somewhere. It consists of three parts:
(i) vector<T> of digits.
(ii) A bool specifying sign of the number.
(iii) A decimal point displacement. (exponent) - Real Algorithm: This data-type is used when a number is represented through a function, it is specially designed for irrational numbers which can be represented by some mathematical function. It consists of three elements:
(i) A function or a lambda function that returns n-th digit of the number.
(ii) A bool specifying sign of the number.
(iii) A decimal point displacement. (exponent) - Real Operation: A number that represents an operation between two real numbers. It consists of:
(i) Two boost::real pointers to its operands.
(ii) An enumerate indicating the type of operation. (+,-,* or /)
These three representations were used to represent a real number in this library. Every operation between two real numbers was represented by a real_operation number, and a real_operation number was a tree pointing to its operands and storing the type of operations it is associated with. Laouen implemented all these types and led the foundation of this library, detailed information on working and implementation of these types can be found at the blog written by Laouen. After that, in GSoC’19, a lot of improvements were done and a basic division algorithm was implemented. A new base was tried and implemented but it was not stable at the time. A back-end class exact_number was created which represents a finite real number, it is used as a back-end number to represent bounds of the interval which represents our real number according to precision and later real_explicit was also modified to store a number in an exact_number, so the members real_explicit class were not a vector, exponent, and bool, it was simply an exact_number which stores a vector, exponent, and bool to represent a finite real number.
After all work done on this library for two years, some improvements were needed to make it ready for the final review. Division and multiplication algorithms were needed to be improved, stable implementation of a larger base was needed, and several more mathematical functionalities were needed to make it more promising. My proposal had the following ideas:
- Implementing a base of 2^n, where n is the number of bits of the variable whose vector will be used to represent the number.
- A separate system for integer and rational number, which will not have a tree structure for operation between an integer and a rational number.
- Implement an integer power method for real numbers.
- Factorial method for integer numbers. (Useful for Taylor expansions).
- Implementing some mathematical functions using Taylor Expansions. (exponential, logarithm, and trigonometric functions.)
- Implementing the Real^Real method, using logarithm and exponential functions.
- Implementing an algorithm for Pi. (Chudnovsky Algorithm.)
- Adding some constants to the library.
Two students were selected to work on this library, me and Kishan Shukla. Although our ideas were independent, there was a huge overlap of ideas in our proposals. So, a lot of work was divided between us. After the division of work, I was left with the following tasks.
- Deciding and implementing a stable and efficient base system. This task was for both of us, me and Kishan, and we planned to complete it during the community bonding period.
- Implement a digit by digit base conversion method for real::algorithm number. (This was implemented during the community bonding period.)
- Introduce integers and rationals as special cases of real number for better optimization possibilities.
- Implement a factorial method for integers.
- Implement Taylor expansions.
- Implement real^real function.
- Add proposed constants.
During community boding, we decided and implemented a stable and efficient base. We agreed on implementing a base of (std::numerical_limits<T>::max()/4)*2 + 1, where T is the type of integral data used to store numbers. (Like int, long, etc).
Then I implemented a digit by digit base conversion method for real::algorithm numbers. Currently, that function is stored as a convenience function, which can be to convert the base of real::algorithm number.
GSoC Phase 1:
Main Goals:
- Documentation for digit by digit base conversion method.
- Implement the proposed system for integer and rational numbers.
- Make tests and documentation for the above system.
Digit by digit base conversion:
This function was added to help when function (for real_operation number) available to the user gives a number that is represented in a different base than internal base representation. For example, if the base of the number that is represented by function is decimal then this function can be used to create another wrapper which will do digit by digit base conversion and return that number in required base.
The function was implemented by representing a digit addition as a multiplication and addition operation. Just assume there is a number N represented in base B1 as N1 and in base B2 as N2(the value of both N1 and N2 are the same, but they are represented in different base). Now, if we append a digit to this number N1 (in base B1), and now we want to produce an equivalent change in N2 (which is stored in base B2), and we are not allowed to do the whole base conversion, we need to do that using N2, both bases B1 and B2, and value of digit that was added in base B1.
N1 = a1a2a3…..an (a1 -> first digit, a2->second digit,…..)
N2 = b1b2b3…..bn(b1-> first digit, a2->second digit,……)
here N1 = N2, but the base of both numbers are different, B1 and B2 respectively. Now we append a digit a(n+1) in N1, what change should we make in N2, to make N2 = N1?
Now, appending digit a(n+1) is equivalent to :
N1 = N1*B1 + a(n+1).
So, we will make the same change in N2:
N2 = N2*B1 + a(n+1)
And that will give us the result of adding a digit in a different base.
Link to initial commit of convenience function: https://github.com/BoostGSoC20/Real/commit/961257de476e5de99fe195c5c2f5f17807b43aa6 (Later the name of the function was modified to append_digit from add_digit).
Integer and rational numbers:
Every operation we do in this library is represented by real_algorithm, which points to numbers on which operation is done and the type of operation which needs to be done. Internally a new copy of operands will be generated and operation will point to those copies. Now, sometimes for small and repetitive operations, the user may not want to create extra copies and a complex tree of operations. Now, a special system for integer and rational number was proposed to avoid these operation trees, any operation involving an integer type and a rational type number will create an integer or rational type number, all other operations like real_explicit + real_rational will create an operation tree. So, they provide some optimization scopes for users.
A class “integer_number” was added to store integer type numbers and another class “real_rational” which stores two “integer_number” objects to represent a rational number in form of “a/b”, where a and b are integers, was added. Initially, integers type numbers were stored by “integer_number” class, and rational numbers were stored by “real_rational” class, later it was decided to use integer_number only in the back-end, to support real_rational class. Then integer numbers(defined by the user) were also represented by real_rational class(in form of a/1). This shifting was done to provide better maintenance of code.
Later, tests for these classes were designed.
GSoC Phase 2:
Main Goals:
- Implement a factorial method for integers.
- Implement Taylor expansions of exponent, logarithm, sin, cos, tan, sec, co-sec, and cot functions.
Factorial:
The main reason behind implementing factorial function was to implement Taylor expansion. We need factorials for Taylor expansions.
The factorial of a number is calculated through repeated multiplication. But after some discussion, this function was not merged into the library. Because, in Taylor expansion, calculating factorial by using factorial of the previous number was more efficient. So, using this function was inefficient, because it calculates every factorial from start. And only use of this function was to help in Taylor expansion, and it was not used in those calculations. So, this function was not added.
The factorial function will calculate factorial of N like this:
n! = 1×2×3×…..×n
but if we have (n-1)! then:
n! = n × (n-1)!
In the calculation of Taylor expansions, we need to calculate (n-1)! or (n-2)! before n! so keeping those factorials and just using them again is better instead of calculating it again from start.
Taylor Expansion:
A Taylor expansion is an infinite converging series, which can represent a continuous and differentiable function, more information can be found on this Wikipedia page.
The pseudo-code of the exponential function:
/**
* EXPONENT FUNCTION USING TAYLOR EXPANSION
* @brief: calculates exponent of a exact_number using Taylor expansion
* @param: num: the number. whose exponent is to be found
* @param: max_error_exponent: Absolute Error in the result should be < 1*base^(-max_error_exponent)* @author: Vikram Singh Chundawat
**/// I am using double here to represent our numbers, but in library exact_number was used, and actual function is a bit different than this. This code is to just explain working of exp function.double exponent(double num, size_t max_error_exponent){double result = 1;
int term_number = 1;
long long factorial = 1;
int cur_term = 0;
double max_error = power((double)1, -max_error_exponent);
double x_pow = 1;
do{
result += cur_term;
factorial *= term_number;
++term_number;
x_pow *= num;
cur_term = x_pow/factorial;
}while(abs(cur_term) > max_error);return result;
}
Four Taylor expansion functions were implemented: exp, ln, sin, and cos. Also, I implemented other trigonometric functions that could be calculated using combinations of these.
Minima and Maxima
The above function is to calculate the exp of a single exact number. But in intervals, we need to check for points of local minima and maxima to show correct results.
Logarithm and exponent do not have local maxima-minima, and they are always increasing function, so exp([a, b]) = [exp(a), exp(b)] and ln([a, b]) = [ln(a), ln(b)].But trigonometric functions are different. Look at the graph of sin(x):
If the input interval [0,π], then it means, our result can be anywhere between 0 and 1. So, the output interval should be [0,1]. But if we directly calculate values of lower and upper bound of input interval, then sin(0) = 0 and sin(π) = 0, then output will be [0,0]. In such cases, where there exists minima or maxima of function between upper and lower bound of input, direct values of the function at upper and lower bound of input interval will give us output which will not contain the correct answer, it will give us wrong output interval.
Trigonometric functions have intervals of 2π, so there can exist a point of minima and maxima between upper and lower bound of input interval, so we need to check for those points.
One way would be to check for multiples of π and π/2 in input interval, but Pi is an irrational number, and checking for its multiples is too inefficient. The first derivative of the function is zero on point of minima and maxima, and it changes its sign at those points. So, points of minima and maxima can be checked using derivatives of functions, we can check whether the sign of derivative is changing or not in input interval. And, we can calculate the derivative of a trigonometric function is also a trigonometric function and we can calculate all trigonometric functions using a combination of sin and cos. So, a function that can calculate both sin and cos in one go was created, it is more efficient than the separate calculation of both functions. Here is the pseudo-code of sin_cos function:
/**
* SINE AND COSINE FUNCTION USING TAYLOR EXPANSION
* @brief: calculates cos(x) and sin(x) of a exact_number using Taylor expansion
* @param: x: the number representing angle in radian
* @param: max_error_exponent: Absolute Error in the result should be < 1*base^(-max_error_exponent)
* @param: upper: if true: error lies in [0, +epsilon]
* else: error lies in [-epsilon, 0], here epsilon = 1*base^(-max_error_exponent)
* @return: a tuple containing sin(x) and cos(x)
* @author: Vikram Singh Chundawat
**/// I am using double here to represent our numbers, but in library exact_number was used, and actual function is a bit different than this. This code is to just explain working of exp function.std::tuple<double, double> sin_cos(double x, size_t max_error_exponent){double sin_result("0");
double cos_result("0");
double cur_sin_term = x;
double cur_cos_term("1");
double cur_power = x;
long long factorial("1");
long long factorial_number("1");
unsigned int term_number_int = 0;
double max_error = pow((double)1, -max_error_exponent);
do{if(term_number_int % 2 == 0){
sin_result += cur_sin_term;
cos_result += cur_cos_term;
}
else{
sin_result -= cur_sin_term;
cos_result -= cur_cos_term;
}
++term_number_int;
factorial_number = factorial_number + literals::one_exact<T>;
factorial *= factorial_number;
cur_power *= x;
cur_cos_term = cur_power/factorial;++factorial_number;
factorial *= factorial_number;
cur_power *= x;
cur_sin_term = cur_power;
cur_sin_term/factorial
}while( (abs(cur_cos_term) > max_error) || (cur_sin_term) > max_error );return std::make_tuple(sin_result, cos_result);
}
Using this function, we can calculate all trigonometric functions, so derivatives can be easily calculated and we can check for points of minima and maxima, and can give results accordingly.
Still, simply using derivatives for maxima-minima points may lead us to the wrong answer. Just consider the input [(-π/2) + e, (3π/2) +e](e is very small value) and we want to calculate sin(x) for this interval. First, we will check for sign of derivative, which is cos(x), the sign of cos(x) will be same at both upper and lower bounds of input interval, so our system will say no maxima-minima points so the result we will be something like [-0.99, -0.98]. If we look carefully at the graph of sin(x), there are both maxima and minima between the upper and lower bound of this input, so our system should give [-1, 1].
The reason for such checking is that a minima-maxima point always alters the sign of derivative. So, the two maxima-minima points will alter the sign of the system two times, so in the end, there will be no change in sign of derivative.
Initially, the system was implemented just like previously proposed, we just checked sign of derivative, later in the third phase we realized this will be wrong and we agreed to make some proper and detailed checks. First, we classified our input between these three categories.
- Upper bound — lower bound < 4.
- 4 ≤ Upper bound — Lower Bound <8.
- Upper bound — Lower bound ≥8.
Someone may wonder why didn’t we use multiples of π for such checking, the reason was the inefficiency of the algorithm used to generate π. This simple check will make this function too slow so we used 4 instead. I am writing here the pseudo-code of the algorithm used to calculate sin(x), which will make things more clear and easy to understand and I am also explaining everything in comments of code. Please go through it:
auto [sin_lower, cos_lower] = sin_cos(lower_bound);
auto [sin_upper, cos_upper] = sin_cos(upper_bound);
/**
First we ensure that our input interval is greater than 2π or not. We can either check that by comparing the difference of upper and lower bound with 2π or a number greater than 2π. We will check whether the difference is greater than 8 or not. If the difference is greater than 8 then we are sure that difference is greater than 2π. We are not using π here for checking because algorithm for π is very complex and checking with that number will be very inefficient. The condition of difference between 2π and 8 will be checked in next case. Here we will only check whether our number is greater than 8 or not. If it is, then we will give hard-coded output of [-1, 1] because the result can be anything between -1 to 1.
**/
if(upper_bound - lower_bound >= 8){
result.lower_bound = -1;
result.upper_bound = 1;
}
/**
Now we will check whether the difference is greater than 4 or not.
4 is greater than π, so we are sure that there is at least one minima/maxima. But there can exist both minima and maxima, as input can be anywhere between [4,8). So, we will check sign of derivative of sin, which is cos. If there is one or three minima/maxima, sign of derivative will change once, so we will different signs of derivative on upper and lower bound. If there are both minima and maxima, sign of derivative will change twice, so at the end, sign of derivative in both upper and lower bound will remain same.
**/
else if(upper_bound - lower_bound >= 4){
/**
If sign of derivative, which is cos(x), if it is same for both upper and lower bound. Then we will return hard-coded output [-1, 1]. Because we are sure that there is at least one minima-maxima, if sign does not changes then it simply means there were both maxima and minima.
**/
if(cos_upper.positive == cos_lower.positive){
result.lower_bound = -1;
result.upper_bound = 1;
}
/**
Now, the sign was changed, then there is either one minima/maxima or three maxima minima. There can exist either one minima/maxima or three points of maxima. In case of three points of maxima, and difference between inputs less than 8, both ends of sin will have same sign and sin(x) of mid point of upper and lower bound will have opposite sign from end points. In that case, we will give hard coded output [-1, 1].
**/
else{
auto mid = (upper_bound + lower_bound) / 2;
if(sin_lower.positive == sin_upper.positive && sine(mid).positive == sin_upper.positive){
result.lower_bound = -1;
result.upper_bound = 1;
}
/**
We will check sign of lower bound of derivative, if it is positive, then initially function was increasing then it would have gone up to 1 and then started decresing. So, the output will be [min(sin(upper_bound, sin(lower_bounf))), 1].
**/
else if(cos_lower.positive){
result.lower_bound = min(sin_upper, sin_lower);
result.upper_bound = 1;
}/**
Now if the sign of derivaye is negative, then initially function was decreasing. So, the function would have got to -1 and then again started incresing. So, the output should be [-1, max(sin_lower, cos_lower)]
**/
else{
result.lower_bound = -1;
result.upper_bound = max(sin_upper, sin_lower);
}
}
}
/**
Now, the final case, the difference is less than 4. There are two possibilities, no maxima/minima or one maxima/minima. If sign of derivative changes, then one maxima/minima, if it doesn't, then no maxima/minima.
**/
else{
/**
If sign does not changes, then no maxima/minima or both maxima and minima. We will check sign of derivative at mid of interval, if it is not same as that of lower and upper bounds. Then we have both maxima and minima.
**/
if(cos_upper.positive == cos_lower.positive){
/* if sign of derivative at mid point is not same as at end points then both minima and maxima are there
*/
auto mid = (upper_bound + lower_bound)/2;
if(cos(mid).positive != cos_lower.positive){
result.lower_bound = -1;
result.upper_bound = 1;
}
/* If cos(x), which is derivative of sin(x), is positive, then function was increasing in that interval. So, the output will be [sin_lower, sin_upper].
*/
else if(cos_lower.positive){
result.lower_bound = sin_lower;
result.upper_bound = sin_upper;
}
/* If it is negative, then function was decreasing in that interval. SO, the output will be [sin_upper, sin_lower].
*/
else{
result.lower_bound = sin_upper;
result.upper_bound = sin_lower;
}
}
/**
If sign changes, then one maxima/minima. So, we will check sign of derivative at lower bound of input. If it is positive, then function was increasing initially. So, the output will be [min(sin_lower, sin_upper), 1]. If it is negative, then function was decreasing initially. So, the output will be [-1, max(sin_upper, sin_lower)].
**/
else{
if(cos_lower.positive){
result.lower_bound = min(sin_lower, sin_upper);
result.upper_bound = 1;
}
else{
result.lower_bound = -1;
result.upper_bound = max(sin_lower, sin_upper);
}
}
}
Theoretically the minima and maxima for all trigonometric functions other than sin and cos are-inf and +inf, we did not have any representation for infinite values and also, returning an interval like [-inf, +inf], does not make sense.
So instead of returning any results, we modified the code to keep iterating for more accurate interval until we get an interval which does not contain any minima or maxima point. If maximum precision is reached and we still don’t get such interval, we throw an error. Pseudo Code for calculation of tan(x) is shown below:
// we will keep on iterating until we get our interval in domain of tan(x)
while(true)
{
auto [sin_lower_tmp, cos_lower_tmp] = sin_cos(input.lower_bound);
auto [sin_upper_tmp, cos_upper_tmp] = sin_cos(input.upper_bound);
bool iterate_again;
if(input.upper_bound - input.lower_bound >= literals::4){
iterate_again = true;
}
else{
/**
Now if difference between lower and upper bounds of interval is less than 4, then there can exist 0,1 or 2 minima/maxima points. First we will check whether the sign of cos(x) from lower to upper bound is changed or not, if it is, then we have one point of maxima/minima. Then we will iterate for better input.
**/
if(cos_upper_tmp.positive != cos_lower.positive || cos_lower_tmp == 0 || cos_upper_tmp == 0){
iterate_again = true;
}
/**
Now, if sign is not changed, then we will check if sign of mid point of derivative is same as end points. If it is, then no points of minima/maxima, no need to iterate further.
**/
else{
auto mid = (input.lower_bound + input.upper_bound ) / 2;
if(cos(mid).positive != cos_lower_tmp.positive){
iterate_again = true;
}
else{
iterate_again = false;
}
}
}
// if we have point of maxima of minima in our input interval
if(iterate_again){
// updating the boundaries of lhs
if(_precision >= input.maximum_precision()){
throw max_precision_for_trigonometric_function_error();
}
input.iterate_n_times(1);
++_precision;
}
else{
sin_lower = sin_lower_tmp;
sin_upper = sin_upper_tmp;
cos_lower = cos_lower_tmp;
cos_upper = cos_upper_tmp;
break;
}
}result.lower_bound = sin_lower / cos_lower;
result.upper_bound = sin_upper / cos_upper;
break;
The logarithm function is not defined for non-positive numbers and sometimes our input can have an interval whose lower bound is negative or zero and upper bound is positive. Then we will iterate for more accurate intervals and show results if the whole interval turns out to be positive else throw an error.
GSoC Phase 3:
Main Goals:
- Implement real^real operation(power operation), using exp and log function.
- Define some commonly used scientific constants.
- Implement tests.
- Document the work.
Real^Real:
The general syntax of a^b can be evaluated using a^b = exp(b * log(a)). We will use the same idea to find non-integral powers of real numbers, integral powers will be calculated using the integer power method created by Kishan Shukla.
Non-integral Powers of negative numbers:
Non-integral powers of negative numbers is a complex number. Its representation is of form (a+ib) form. More information can be found here. Currently, our library does not support imaginary numbers, so we throw an error for such inputs.
Here is the code for power function:
power(real_num, power){
// checking whether the number is integer or not
try{
result = real(real_operation<T>(real_num, power, OPERATION::INTEGER_POWER));
return result;
}
catch(const negative_integers_not_supported& e1){
power *= -1;
result = real(real_operation<T>(real_num, power, OPERATION::INTEGER_POWER));
result = real(real_operation<T>(one, result, OPERATION::DIVISION));
return result;
}
catch(const non_integral_exponent_exception& e2){
try{
/**
* result = exp(exponent * log(real_num))
* Warning: the result of non-integral power of a negative number is a complex number, and we do not support complex numbers. So, that will generate error.
* Note: The error will be generated by logarithm function.
**/
/**
* Now, if number is negative, then logarithm function will check out and throw error
**/
result = real(real_operation<T>(real_num, zero, OPERATION::LOGARITHM));
result = real(real_operation<T>(result, power, OPERATION::MULTIPLICATION));
result = real(real_operation<T>(result, zero, OPERATION::EXPONENT));
return result;
}
catch(const logarithm_not_defined_for_non_positive_number& e3){
throw non_integral_power_of_negative_number();
}
}
}Scientific constants:
Some common constants like mass and charge of subatomic particles were added. The values were taken from the most accurate measurements available.
Square root function:
Although this function was not in my proposal, but I was left with some time in the end. So I added this function to the library. This function simply calculates x^(1/2) using logarithm and exponent function.
After that, all remaining tests were made and some minor bugs which were detected at the time of testing were fixed. After that, this report was written along with some documentation in code.
Doxygen documentation will be automatically generated from comments written in code.
This entire project was a very good journey and we learned a lot during this period, and also understood a lot of standard practice that is followed by programmers when dealing with big projects. Special thanks to mentors Damian Vicino and Laouen Belloli. They constantly gave us feedback on our work and taught us a lot.
Future Work
After finishing all work that was necessary to make this library ready for review, we need to do some optimizations. This library always makes iterations from start for every calculation, this leads us to many repeated calculations. One can optimize this by storing the number at its last calculated precision. One more important optimization, also proposed by Laouen in his report, replace Tree Representation of the library to Directed Acyclic Graph(DAG).
