A Problem of Accuracy of Computer Calculations

The paper presented the results of the research related to the analysis of the reliability of computer calculations. Relevant examples of incorrect program operation were demonstrated: both quite simple and much less obvious, such as S. Rump's example. In addition to mathematical explanations, authors focused on purely software capabilities for controlling the accuracy of complex calculations. For this purpose, examples of effective use of the functionality of the decimal and fraction modules in Python 3.x were given.

The problem of reliability (accuracy) of computer calculations is one of the fundamental problems of computer science [1, 2]. Physical limitations on the amount of memory allocated by computer technology for number processing impose fundamental limitations on the possibilities and logic of organizing computer calculations, which also requires the involvement of specific mathematical research methods.

At the same time, the problem is not only the possibility of an incorrect result of theoretical research. As practice shows, the consequences of incorrect computer calculations can be catastrophic. There is a famous case: on February 25, 1991, during the Gulf War I (Operation Desert Storm), a “Patriot” air defence system failed to intercept an enemy missile, and the projectile hit a barracks of US soldiers, killing 28 people. An official investigation [3] showed that 24-bit interceptor processors make a time conversion error of 0,013 seconds every hour. “Patriot” was not restarted for more than 100 hours, and as a result of the accumulation of such seemingly minor errors, during this time there was an error in calculating the position of the missile at 600 meters.

Almost any academic textbook about calculation methods – an integral discipline of the training course for both physical and mathematical specialists and computer science specialists – usually contains only mathematical information about errors and the accuracy limits of calculations. The applied side of the issue, i.e. the computer implementation of calculations, remains out of consideration. Technical disciplines, such as “Programming”, usually provide only the basic concepts of converting numbers to a binary code and back, and provide only general information about the existing standards for representing floating-point numbers.

Ignoring the fundamentals of binary arithmetic and the principles of the IEEE-754 standard [4] can significantly affect the correctness of both theoretical research (if software or computer equipment was used for modelling) and simply the quality of the product when developing application programs.

Programs are implemented in C++, Java, or Python programming languages, depending on when it is more appropriate: in C++ and Java, we can visually experiment by changing the types of floating-point numbers, float and double; in Python, by default, “long” arithmetic is implemented, that is, a set of algorithms for bitwise work with numbers of arbitrary length.

Consider the following example of a feature of computer calculations. Let's check with the help of the C++ programming language, whether the distributive law is fulfilled in computer calculations. Two simple program codes for this can be seen in Figure 1 below.

Figure 1

Program codes in the C++ programming language for checking of the fulfilment the distributive law: (a) for numbers 0,987, 0,278, and 1,268 in double format; (b) for numbers 0,9871, 0,2781, and 1,268 in float format

As expected, for code in Figure 1 (a) we get a result of $1$ (i.e. $\mathrm{true}$).

Let's slightly change the values of the numbers a and b, increasing them by 0,0001 (that is, we will change the numbers a and b by $\mathrm{0,01}\%$ and $\mathrm{0,04}\%$, respectively). As a result, we will get a result of $0$ ($false$).

If we output to the console the corresponding values for the case $a=\mathrm{0,9871}$ and $b=\mathrm{0,2781}$ with an accuracy of 20 decimal places, the program will issue:

$a\cdot c+b\cdot c\to 1,60427359999999996631,$

$\left(a+b\right)\cdot c\to 1,60427360000000018836.$

So we get different values, which are insignificant but differ from the real value of $\mathrm{1,6042736}$ (the absolute error is $3,4\cdot {10}^{-17}$ and $1,9\cdot {10}^{-16}$).

Let's change the types of variables in the program to float (Figure 1 (b)), and as a result, we will already get the result $true$. That is, it is the reduction in the accuracy of the calculations that returns us to the correct result from a mathematical point of view, which at first glance looks somewhat paradoxical.

The above complexities of computer calculations arise from the way real numbers are represented in computer memory, regulated by IEEE-754 standards.

The description of the IEEE-754 standard goes beyond the scope of this work, so we will limit ourselves to the demonstration of the representation of numbers from the example described above in float and double formats.

Float format: 32 bits are allocated to the number: the first bit is the sign bit (if the sign bit is 0, the number is non-negative (positive or zero); if the sign bit is 1 then the number is negative), the next 8 bits are allocated to the exponent, and the last 23 bits are allocated to the normalized mantissa. This is shown schematically in Figure 2 below.

Figure 2

Representation of a number in float format in the IEEE-754 standard

Double format (i.e. double-precision format): 64 bits are allocated to the number: the first bit is the sign bit – the bit that indicates the sign of a number; 11 bits are allocated to the exponent, and the last 52 bits are accordingly allocated to the normalized mantissa (Figure 2).

Figure 3

Representation of a number in the double format in the IEEE-754 standard

A simple example of comparing these two formats is comparing the same number in different formats using the Java 8 programming language.

Figure 4

The program for comparing of the number 0,3 presented in different formats

This example compares the number 0,3 represented in float and double formats.

Figure 5

The result of the program from Figure 4

As a result, we can see that the absolute values of two seemingly identical numbers are not equal to each other and have a significant difference when it comes to the accuracy of calculations.

Next, we will consider the example proposed by S. Rump [5] back in 1988. It is necessary to calculate the value of the expression

$f\left(x,y\right)=333,75\cdot y^6+x^2\cdot \left(11\cdot x^2{y}^2-y^6-121\cdot y^4-2\right)+5,5\cdot y^8+\frac{x}{2y}$

at the point $(77617 ; 33096)$.

The task looks quite simple. If we temporarily ignore the last term $\frac{x}{2y}$ (which is not really important), and considering that $y = 33096$ is an even number, it is obvious that the value of the expression $f\left(x,y\right)-\frac{x}{2y}$ in general must be an integer.

To implement the Rump’s example, we will use Python 3

Figure 6

The program implementation of the Rump’s example

The result of executing the program in Figure 6 will be the value $1180591620717411303424,$ although in fact

$f\left(77617,33096\right)=-\frac{54767}{66192}$

The characteristic problem is precisely at the point $(77617;33096).$

Let's explain what is special about this point. Let's separate the positive and negative parts of the expression (as already mentioned: the last addition of $\frac{x}{2y}$ does not really significantly affect the correctness of the calculation result). Let's separate the positive and negative parts of the expression $\left(x,y\right)-\frac{x}{2y}$ , as shown in Figure 7.

Figure 7

The program implementation of the Rump’s example with separated positive and negative parts

We will get the result

$positive=7917112216566288664689761316426849984;$

$negative=7917112216566288664689761316426849986.$

That is, in our example

$positive–negative=-2,$

so final result coincides with the mathematical one. The fact is that one of the main problems when working with floating-point numbers is the subtraction operation, because when subtracting numbers close in value, significant digits may be lost, and the result in this case will be incorrect.

If we change one line in the program by replacing integer division (//) with regular division (/), that is

$positive=\frac{1335\cdot y^6}{4}+x^2\cdot \left(11\cdot x^2{y}^2\right)+\frac{11\cdot y^8}{2}$

so the value of the variable 𝑝𝑜𝑠𝑖𝑡𝑖𝑣𝑒 will be equal to

$positive=7917112216566288629721075632129966080$

and the difference between the values of 𝑝𝑜𝑠𝑖𝑡𝑖𝑣𝑒 and 𝑛𝑒𝑔𝑎𝑡𝑖𝑣𝑒 will be already

$-34968685684296883906.$

If we rewrite the expression for the variable $positive$ in a format more familiar to the Python user

$positive=\frac{1335}{4}\cdot y^6+x^2\cdot \left(11\cdot x^2\cdot y^2\right)+\frac{11}{2}\cdot y^8,$

then the value of 𝑝𝑜𝑠𝑖𝑡𝑖𝑣𝑒 will be equal to

$positive=7917112216566289810312696349541269504$

since different rounding errors are given in signs for different arithmetic operations (as we already noted in the previous section: the order of operations also matters for the accuracy of calculations). The difference in this case between the received positive value and the previous one will be $1180591620717411303424,$ and the value

$positive–negative=1145622935033114419518.$

Correctly change the sequence of performing arithmetic operations and immediately get errors of the order ${10}^{21}$ .

It should be noted that Python has the ability to a certain extent to eliminate some problems arising from calculations with floating-point numbers – for example, using the functionality of the decimal module to control the accuracy of calculations. The official documentation for the decimal module can be found at [6]. Without going into unnecessary details, we note that the decimal module (in particular, the Decimal class) makes it possible to work with floating-point numbers as with ordinary fractions.

With all of the above, for example, only using the capabilities of the Decimal class seems insufficient. Skipping several intermediate stages of the research, the author found a way to get the correct result of software calculations. This method consists in the joint use of the decimal module together with the application of the multiplication reduction scheme (Horner scheme), with the aim of reducing the number of arithmetic operations and increasing the accuracy of calculations. In this case, we have to rewrite the expression $f\left(x,y\right)-\frac{x}{2y}$ in the form

$f\left(x,y\right)=y^2\cdot \left(y^2\cdot \left(y^2\cdot \left(5,5\cdot y^2\left(333,75-x^2\right)\right)-121\cdot x^2\right)+11\cdot x^4\right)-2\cdot x^2.$

We implement the described approach programmatically (Figure 8).

Figure 8

The program implementation of the example using the reduced multiplication scheme

In the result we get expected value $\mathrm{-2}$.

The described examples can be called linear, since the reasons for the occurrence of calculation features become obvious even with a cursory analysis, if we take into account the IEEE-754 standard.

The question of using the results mentioned in the paper, in particular, in the modelling of random and natural processes was considered in recent studies of such processes [7, 8, 9, 10, 11, 12, 13, 14, 15, 16].

Author Contributions: All authors have read and agreed to the published version of the manuscript.

Funding: This research received no external funding.

Acknowledgments: The authors express their heartfelt gratitude to the brave soldiers of the Ukrainian Armed Forces who protect the lives of the authors and their families from Russian bloody murderers since 2014.

Conflicts of Interest: The authors declare no conflict of interest.