Error Compensation for Floating-Point Operations

The work is devoted to rounding errors arising in the calculations of floating-point numbers. Here the following topics will be briefly considered: “Representation of real numbers”, “Methods for finding rounding errors in floating point numbers” and an example of compensation for rounding errors will be given.

In this paper, examples are given in the programming language C.

Representation of real numbers

Consider the representation of a finite real number in the IEEE 754-2008 standard as an expression that is characterized by three elements: S (0 or 1), the mantissa M, and the order E :

v = -1 ^S * b ^{(E - BIAS)} * M

• base b (the standard defines three binary formats and two decimal);
• the length of the mantissa p , which determines the accuracy of the representation of the number, is measured in numbers (binary or decimal);
• maximum order value emax ;
• minimum order value emin ; the standard requires that the condition:
  

   emin = 1-emax 

• BIAS (offset); order is written in the so-called. biased form, i.e. the
  real value of the exponent is equal to E-BIAS .

The table below shows the parameters of standard floating-point number formats. Here: w is the width of the bit field to represent the order, t is the width of the bit field to represent the mantissa, k is the total width of the bit string.



The following table shows the ranges and accuracy of the standard 32 and 64-bit floating point real-life formats.



Here "Precision, Epsilon" is the smallest number for which the expression is true:

1 + EPSILON != 1

This Epsilon value characterizes the relative accuracy of addition and subtraction operations: if the value added to x or subtracted from x is less than epsilon * x , then the result will remain equal to x . In practice, in some cases, when using quantities approaching epsilon * x in additive operations, round-off errors of a smaller term begin to affect . Such situations will be discussed in this paper.

Consider the representation of a real number for the data type float (Binary32).



In this example:

• Number V = 0.15625 ₁₀
• Sign S = 0, i.e. +
• Order (E - BIAS) = 01111100 ₂ - 01111111 ₂ = 124 ₁₀ - 127 ₁₀ = -3 ₁₀
• Mantissa M = 1.01000000000000000000000 ₂

Thus, the number V = 1.01 ₂ * 2 ₁₀ ^{-3 ₁₀} = 101 ₂ * 2 ₁₀ ^{-5 ₁₀} = 5 ₁₀ * 2 ₁₀ ^{-5 ₁₀} = 0.15625 ₁₀

Note that when performing operations with real numbers, often the result will not fit in N-bit matises , i.e. rounding will occur. Rounding in one computational operation does not exceed the order of EPSILON / 2 , but when we need to do a lot of operations, to improve the accuracy of calculating the result, we need to learn how to find exactly how the result of each particular operation is rounded.

For more information about rounding and violation of axiomatics, see the file makarov_float.pdf (link to the material below).

Ways to find rounding errors for floating point numbers

This problem was investigated by many experts, the most famous of them are: David Goldberg, William Kachen, Jonathan Richard Shevchuk.

Below we consider the algorithms for finding rounding errors given in the work of Shevchuk, using two functions as an example:

• TwoSum - the function of finding the rounding error when adding.
• TwoProduct - the function of finding the rounding error when multiplying.

For a correct understanding of these algorithms, we will rely on the theorems considered in Shevchuk's work. We give the theorems without proof.

Saying that the number a is a p- bit number, it means that the length of the matisse of a is represented by p bits.

Twosum

Theorem: let the numbers a and b be p- bit floating point numbers, where p > = 3, then following this algorithm we will get 2 numbers: x and y , for which the condition is fulfilled: a + b = x + y . Moreover, x is an approximation of the sum of a and b , and y is the rounding error of the calculation of the number x .

Twoproduct

In fact, the algorithm for finding rounding errors when multiplying two real numbers consists of 2 functions: Split - an auxiliary function and the TwoProduct function, where we find the error.

Consider the algorithm of the Split function.

Split

Theorem: the number a is a p- bit floating point number, where p > = 3. Choose a break point s , where p / 2 <= s <= p -1. Following the algorithm, we get a ( p - s ) -bit number - the number a_hi and ( s -1) -bit number - a_lo , where | a_hi | > = | a_low | and a = a_hi + a_low .



Now let's move on to the analysis of the TwoProduct function.

Twoproduct

Theorem: let the numbers a and b be p- bit floating point numbers, where p > = 6. Then, performing this algorithm, we will get 2 numbers: x and y , for which the condition is satisfied: a * b = x + y . Moreover, x is an approximation of the product of the numbers a and b , and the number y is the rounding error of the calculation of the number x .



It can be seen that the final result is represented by a pair of N-bit real numbers: result + error. Moreover, the second should have an order of not more than EPSILON with respect to the first.

An example of the compensation of rounding errors in floating point numbers

Now we turn to practical calculations, which are based on Shevchuk's algorithms. We find rounding errors in the addition and multiplication of floating point numbers and analyze how the error accumulates.

Sum rounding error

Here is an example of a simple summation program:

#include 
#include 
int main() {
    float val = 2.7892;
    printf("%0.7g \n", val);
    val = val/10000000000.0;
    float result = 0.0;
    for (long long i = 0; i < 10000000000; i++) {
        result += val;
    }
    printf("%0.7g \n", result);
    return 0;
}

As a result of the program, we should get two identical numbers: 2.7892 and 2.7892. But the console output was: 2.7892 and 0.0078125. This shows that the error accumulated is very large.

Now let's try to do the same, but using the Shevchuk algorithm, we will accumulate the error in a separate variable, and then we will compensate the result by adding the error to the main variable of the sum.

#include 
#include 
float TwoSum(float a, float b, float& error) {
    float x = a + b;
    float b_virt = x - a;
    float a_virt = x - b_virt;
    float b_roundoff = b - b_virt;
    float a_roudnoff = a - a_virt;
    float y = a_roudnoff + b_roundoff;
    error += y;
    return x;
}
int main() {
    float val = 2.7892;
    printf("%0.7g \n", val);
    val = val/10000000000.0;
    float result = 0.0;
    float error = 0.0;
    for (long long i = 0; i < 10000000000; i++) {
        result = TwoSum(result, val, error);
    }
    result += error;
    printf("%0.7g \n", result);
    return 0;
}

As a result, we get 2 numbers: 2.7892 and 0.015625. The result has improved, but the error still makes itself felt. In this example, the error that occurs in the addition operation was not taken into account in the TwoSum () function:

error += y;

We will compensate for the result at each iteration of the cycle and overwrite the error in a variable that accumulates the rounding error in the summation operation. To do this, we modify the TwoSum () function: add an isNull variable of type bool which indicates whether we should accumulate the error or whether it should be rewritten.

As a result, result will be represented by 2 variables: result is the main variable, error1 is the error of the operation result + = val .

The code will look like this:

#include 
#include 
float TwoSum(float a, float b, float& error, bool isNull) {
    float x = a + b;
    float b_virt = x - a;
    float a_virt = x - b_virt;
    float b_roundoff = b - b_virt;
    float a_roudnoff = a - a_virt;
    float y = a_roudnoff + b_roundoff;
    if (isNull) {
        error = y;
    } else {
        error += y;
    }
    return x;
}
int main() {
    float val = 2.7892;
    printf("%0.7g \n", val);
    val = val/10000000000.0;
    float result = 0.0;
    float error1 = 0.0;
    for (long long i = 0; i < 10000000000; i++) {
        result = TwoSum(result, val, error1, false);
        result = TwoSum(error1, result, error1, true);
    }
    printf("%0.7g \n", result);
    return 0;
}

The program will output the numbers: 2,7892 and 2,789195.

Note that here the rounding error arising in the multiplication operation was not taken into account:

val = val*(1/10000000000.0);

We take this error into account by adding the functions of accounting for multiplication errors, which were developed by D.R. Shevchuk. In this case, the variable val will be represented by two variables:

val_real = val + errorMult

Thus, result will be represented by 3 variables: result is the main variable, error1 is the error of the operation result + = val , error2 is the error of the operation result + = errorMult .

We will also add the error1 and error2 variables , and write the error from this operation into error2 .

As a result, the code:

#include 
#include 
float TwoSum(float a, float b, float& error, bool isNull) {
//isNull отвечает за то, стоит ли нам накопить возникающую погрешность
// или же стоит перепизаписать ее
    float x = a + b;
    float b_virt = x - a;
    float a_virt = x - b_virt;
    float b_roundoff = b - b_virt;
    float a_roudnoff = a - a_virt;
    float y = a_roudnoff + b_roundoff;
    if (isNull) {
        error = y;
    } else {
        error += y;
    }
    return x;
}
void Split(float a, int s, float& a_hi, float& a_lo) {
    float c = (pow(2, s) + 1)*a;
    float a_big = c - a;
    a_hi = c - a_big;
    a_lo = a - a_hi;
}
float TwoProduct(float a, float b, float& err) {
    float x = a*b;
    float a_hi, a_low, b_hi, b_low;
    Split(a, 12, a_hi, a_low);
    Split(b, 12, b_hi, b_low);
    float err1, err2, err3;
    err1 = x - (a_hi*b_hi);
    err2 = err1 - (a_low*b_hi);
    err3 = err2 - (a_hi*b_low);
    err += ((a_low * b_low) - err3);
    return x;
}
int main() {
    float val = 2.7892;
    printf("%0.7g \n", val);
    float errorMult = 0;//погрешность умножения
    val = TwoProduct(val, 1.0/10000000000.0, errorMult);
    float result = 0.0;
    float error1 = 0.0;
    float error2 = 0.0;
    for (long long i = 0; i < 10000000000; i++) {
        result = TwoSum(result, val, error1, false);
        result = TwoSum(result, errorMult, error2, false);
        error1 = TwoSum(error2, error1, error2, true);
        result = TwoSum(error1, result, error1, true);
    }
    printf("%0.7g \n", result);
    return 0;
}

The following numbers were displayed on the console: 2.7892 and 2.789195.

This suggests that the rounding error of the multiplication is small enough to appear even at 10 billion iterations. This result is as close as possible to the original number, if we take into account the errors in the operations of addition and multiplication. To obtain a more accurate result, you can enter additional variables that take into account errors. Say, add a variable that takes into account the error in the operation of accumulating the basic error in the TwoSum () function. Then this error will be of the order EPSILON ² , with respect to the main result (the first error will be of the order EPSILON ).

Multiplication Rounding Error

We calculate the number 1.0012 ¹⁰¹ in the cycle, i.e. we will do the following:

#include 
#include 
int main() {
    float val = 1.0012;
    float result = 1.0012;
    for (long long i = 0; i < 100; i++) {
        result *= val;
    }
    printf("%0.15g \n", result);
    return 0;
}

Note that the exact result, accurate to the fifteenth decimal place, is 1.128768638496750. We will receive: 1.12876391410828. As you can see, the error turned out to be quite large.

We will derive the variable val, converting it to a double data type, and see what is actually written into it:

printf("%0.15g \n", (double)val);

We get the number 1.00119996070862. This suggests that in programming, even the most accurate constant is neither reliable nor constant. Therefore, our real accurate result will be 1.128764164435784, accurate to the fifteenth decimal place.

Now let's try to improve the result obtained earlier. To do this, we introduce compensation for the calculation result by taking into account the rounding error in the multiplication operation. We will also try to add the accumulated error to the result variable at each step.

The code:

#include 
#include 
float TwoSum(float a, float b, float& error, bool isNull) {
//isNull отвечает за то, стоит ли нам накопить возникающую погрешность
// или же стоит перепизаписать ее
    float x = a + b;
    float b_virt = x - a;
    float a_virt = x - b_virt;
    float b_roundoff = b - b_virt;
    float a_roudnoff = a - a_virt;
    float y = a_roudnoff + b_roundoff;
    if (isNull) {
        error = y;
    } else {
        error += y;
    }
    return x;
}
void Split(float a, int s, float& a_hi, float& a_lo) {
    float c = (pow(2, s) + 1)*a;
    float a_big = c - a;
    a_hi = c - a_big;
    a_lo = a - a_hi;
}
float TwoProduct(float a, float b, float& err) {
    float x = a*b;
    float a_hi, a_low, b_hi, b_low;
    Split(a, 12, a_hi, a_low);
    Split(b, 12, b_hi, b_low);
    float err1, err2, err3;
    err1 = x - (a_hi*b_hi);
    err2 = err1 - (a_low*b_hi);
    err3 = err2 - (a_hi*b_low);
    err += ((a_low * b_low) - err3);
    return x;
}
int main() {
    float val = 1.0012;
    float result = 1.0012;
    float errorMain = 0.0;
    for (long long i = 0; i < 100; i++) {
        result = TwoProduct(result, val, errorMain);
        result = TwoSum(errorMain, result, errorMain, true);
    }
    printf("%0.15g \n", result);
    return 0;
}

The program displays the following number: 1.12876415252686. We got an error of 1.0e-008, which is less than EPSILON / 2 for the float data type . Thus, this result can be considered quite good.

Summary

1) In this paper, we considered the representation of floating point numbers in the format of the IEEE 754-2008 standard.
2) A method was shown for finding rounding errors in addition and multiplication for floating point numbers.
3) Simple examples of compensating rounding errors for floating point numbers were considered.

The work was performed by Victor Fadeev.
Advised Makarov A.V.

PS Thanks for the errors found to users:
xeioex , Albom .

References

• The work of Andrei Makarov on the theme "Representation of real numbers." The makarov_float .pdf file shows how rounding errors occur with floating point numbers.
• Wikipedia . Single precision number.
• Adaptive Precision Floating-Point Arithmetic
  and Fast Robust Geometric Predicates
  Jonathan Richard Shewchuk