The Last Digit

Co-authored by Claude · Version 1.0.0, changed in rev 133

Why an answer can be wrong in its last digit, two times an error of that kind grew until it mattered, and the test DKE Python holds its arithmetic to

Ask Python 3.11 for the logarithm of 1,000 to base 10 and it answers 2.9999999999999996. The answer is three, exactly, and three is a number a computer holds exactly. Python's misses it by about 0.0000000000000004, a difference in the sixteenth decimal place, and it is still not three.

This essay is about that last digit: where it comes from, why the people who write the arithmetic usually let it be, the two times an error of its kind grew until it mattered, and the stricter test DKE Python sets itself. Python's answers here were reproduced on CPython 3.11.2 and CPython 3.12.16, both over the GNU C Library (glibc) 2.36 as Debian 12 ships it. Every exact value was computed independently to at least 80 digits and rounded once. Everything the essay says about DKE is something DKE's published documents say.

A tenth that is not a tenth

A Python float is a binary floating-point number of the kind the IEEE 754 standard defines: 53 binary digits and an exponent, about sixteen decimal digits. Most decimal fractions cannot be held exactly that way. A tenth in binary is 0.000110011001100…, with the 0011 repeating forever, so the number Python keeps for 0.1 is the nearest one it can hold. Python's own tutorial prints it in full: 0.1000000000000000055511151231257827021181583404541015625.

Every arithmetic step has an exact result, and that result usually falls between two numbers the machine can hold, so the step has to pick one. Since its first edition in 1985, IEEE 754 has required addition, subtraction, multiplication, division and the square root to be correctly rounded. The answer is the exact result of the operation on the numbers given, rounded once to the nearest number that can be held. When the exact result lies exactly halfway, the answer is the neighbour whose last binary digit is even.

That is why 0.1 + 0.2 printing 0.30000000000000004 is not an error. The two numbers added are not quite a tenth and two tenths, and their exact sum rounds to the number that prints that way. The arithmetic did what it promised with the numbers it was given. The surprise was in the numbers.

Where the promise stops

The promise covers those few operations. It does not cover exp, log, sin or cos. The 2008 revision of IEEE 754 lists them and many more in its Clause 9 as recommended operations, which must round correctly where a language provides them. None of them is required in order to conform to the standard.

There is a reason for the gap. To round an answer correctly you have to know which side of a halfway point the exact value lies on, and an exact value can lie very close to one. William Kahan, the key designer of IEEE 754, named the difficulty the Table-Maker's Dilemma: "Nobody knows how much it would cost to compute yʷ correctly rounded for every two floating-point arguments at which it does not over/underflow. Instead, reputable math libraries compute elementary transcendental functions mostly within slightly more than half an ulp and almost always well within one ulp." An ulp is a unit of the last place: the gap between a number and the next one the machine can hold.

The GNU C Library supplies these functions to most programs on Linux, and it says the same of itself. In its manual for version 2.36: "the GNU C Library does not aim for correctly rounded results for functions in the math library". Its goal is a result "within a few ulp" of the mathematically correct value. CPython's documentation says that its math module "consists mostly of thin wrappers around the platform C math library functions". So on Debian 12, Python's cosine is the C library's:

import math
print(math.cos(-55.223308363883085))
print(math.log10(1.4142135623730951))

0.24298017990326443
0.1505149978319906

The exact cosine at that argument, rounded once, is 0.24298017990326445, and the exact logarithm is 0.15051499783199063. Each answer is one unit of the last place away, which is within what the library says it aims for. Measured over 200,000 arguments each on that library, exp, log, sin, cos, atan2, hypot and pow missed the correctly rounded answer between about one time in 600 and one time in 5,000, always by one unit of the last place.

Two roundings where one would do

Some last digits are lost not inside a function but between two of them. CPython's documentation says how math.log with a base is computed: "calculated as log(x)/log(base)". Both logarithms are rounded, and then their quotient is rounded again:

import math
print(math.log(1000, 10), math.log(125, 5), math.log(243, 3))

2.9999999999999996 3.0000000000000004 4.999999999999999

The exact answers, 3, 3 and 5, can all be held exactly, and none of the three printed is one of them. math.degrees and math.radians do the same with a constant. They multiply by 180/π or π/180 already rounded to a float, and round the product again.

Complex numbers compose more steps. The imaginary part of cmath.exp(complex(1.3, 2.7)) is e to the power 1.3 times the sine of 2.7. Python computes each factor and rounds it, then multiplies and rounds again. It answers 1.5681835703495035, where the exact value rounds to 1.5681835703495033.

One rounding after another

A sum rounds once for every addition. The Python tutorial's own example, and one more:

print(sum([0.1] * 10), sum([1e16, 1.0, -1e16]))

0.9999999999999999 0.0

Ten roundings leave the first a unit short of one. The second loses the 1.0 entirely: next to 10¹⁶ the numbers a float can hold are 2 apart, so 10¹⁶ + 1 rounds back to 10¹⁶ before the subtraction can recover it.

Python offers a remedy in math.fsum, which its documentation calls "an accurate floating point sum", kept by "tracking multiple intermediate partial sums". Python 3.12 changed sum itself, and its documentation says so: "Summation of floats switched to an algorithm that gives higher accuracy on most builds." On CPython 3.12.16 the line above prints 1.0 1.0.

Higher accuracy is what 3.12 promises, and that is not the same as the exact sum rounded once. Add 1.0, 10⁻¹⁶ and 10¹⁶. The exact total is 10,000,000,000,000,001.0000000000000001. The two floats nearest it are 10¹⁶ and 10¹⁶ + 2, so 10¹⁶ + 1 is exactly halfway between them, and the tiny 10⁻¹⁶ puts the total just over the halfway point. Rounded once, the answer is 10¹⁶ + 2, which Python prints as 1.0000000000000002e+16. Python 3.11 and Python 3.12.16 both answer 1e+16. Each loses the tiny one in a rounding along the way, which leaves the total exactly halfway, and a tie goes to 10¹⁶, the neighbour whose last binary digit is even.

When the last digit was not the last

An error of one unit of the last place is one or two parts in ten million billion. On its own it rarely matters. Twice in the 1980s and 1990s, an error of the same kind was made the same way thousands or millions of times, and grew.

Dhahran, 1991

On 25 February 1991, during the Gulf War, a Patriot missile battery at Dhahran in Saudi Arabia failed to track and intercept an incoming Iraqi Scud missile. The Scud struck an American Army barracks and killed 28 soldiers. The United States General Accounting Office reported the cause in February 1992, in Patriot Missile Defense: Software Problem Led to System Failure at Dhahran, Saudi Arabia (GAO/IMTEC-92-26), and Douglas Arnold's account of the failure works through its arithmetic.

The system's clock counted time in tenths of a second, as a whole number. To predict where a Scud would next appear, it turned that count into seconds, and in the report's words, "because of the way the Patriot computer performs its calculations and the fact that its registers are only 24 bits long, the conversion of time from an integer to a real number cannot be any more precise than 24 bits." The tenth it multiplied by was cut short, not rounded. Arnold writes it out as the register kept it, 0.00011001100110011001100 in binary, which is 838,860 parts in 2²³ and falls short of a true tenth by about 0.000000095. Every tenth of a second the clock had counted carried that shortfall, and the battery had been running for over 100 hours. The report tabulates how far the converted time had drifted after each length of continuous running, and the shortfall reproduces every line of that table, to the four places it prints:

from fractions import Fraction
kept = Fraction(838860, 2 ** 23)
short = Fraction(1, 10) - kept
print(float(short))
for hours in (1, 8, 20, 48, 72, 100):
    print(hours, round(float(short * hours * 60 * 60 * 10), 4))

9.5367431640625e-08
1 0.0034
8 0.0275
20 0.0687
48 0.1648
72 0.2472
100 0.3433

A third of a second. In the same table that drift has moved the range gate, the region in which the radar looked for the Scud, by 687 metres, and a note says that after about 20 hours of continuous running the target was already outside it. In the report's words, the inaccuracy was by then "serious enough to cause the system to look in the wrong place for the incoming Scud." Arnold adds that the time calculation had been improved in some parts of the code but not all, so the errors that might have cancelled did not.

Vancouver, 1982 to 1983

In January 1982 the Vancouver Stock Exchange started a new index at 1,000. It was recalculated with every trade, about 3,000 times a day, and after each recalculation the value was truncated to three decimal places and carried into the next. Each cut was less than a thousandth of a point, and every cut went down. The index lost about 25 points a month that the market had not lost. When the error was corrected over the weekend of 25 to 28 November 1983, the index went from its Friday close of 524.811 to 1,098.892. The Wall Street Journal reported the problem on 8 November 1983, and the Toronto Star the correction on 29 November.

Both failures have the same shape. A tiny error, made the same way every time, was carried forward and added up. Cutting a number short always errs in the same direction, so nothing cancels. Rounding to the nearest errs in either direction, and its errors partly cancel: ten roundings left the ten tenths only one unit short.

A stricter test

None of the Python answers above is a bug. The C library says what it aims for, and every miss measured here is within it. CPython says how it computes log with a base and how it sums, and does exactly that. Each answer is within the contract of the project that gave it.

DKE Python, the language of DKE, the Deterministic Knowledge Engine, holds its arithmetic to a stricter test: the exact value, rounded once. Its Specification states the test function by function:

• A sum over numbers holding a real "is the exact sum of the start and every element, rounded once to the nearest real, ties to even, so the order the elements come in cannot change it".
• "Powers, logarithms and moduli are correctly rounded."
• math.log(x, base) is "the logarithm of x to base rounded once", and math.degrees(x) and math.radians(x) are x times 180/π and π/180, "each exact product rounded once".
• cmath.exp, cmath.rect and cmath.sqrt answer "each part of the result, the real and the imaginary, as its exact value rounded once".

The same calls in DKE Python:

import math
import cmath

print(sum([0.1] * 10), sum([1e16, 1.0, -1e16]), sum([1.0, 1e-16, 1e16]))
print(math.log(1000, 10), math.log(125, 5), math.log(243, 3))
print(math.cos(-55.223308363883085), math.log10(1.4142135623730951))
print(math.degrees(-880.7976600675347))
print(cmath.exp(complex(1.3, 2.7)).imag)

1.0 1.0 1.0000000000000002e+16
3.0 3.0 5.0
0.24298017990326445 0.15051499783199063
-50465.9885268683
1.5681835703495033

What the stricter test costs is agreement with Python in the last digit. DKE Python follows Python's own algorithm for each function, and the Specification says where it departs and why: "Seven functions depart from Python's algorithm, because it rounds where the answer need not be." It is plain about the one thing it cannot promise. Where an exact value lies so near the midpoint between two reals that it could round either way, the call "is not proven to round it correctly, and no such case is known."

What it gives is an answer that, in the Specification's words, "depends on its operands alone, and is the same wherever it is computed". Python's cosine is whatever the C library under it computes. DKE Python's is the exact cosine, rounded once.

What rounding once cannot do

Rounding once removes the error a calculation adds. It cannot remove the error already in its numbers:

s = sum([0.1, 0.1, 0.1])
print(s, s == 0.3)

0.30000000000000004 False

The three stored tenths add up, exactly, to a number precisely halfway between the real written 0.3 and the next one up. A tie goes to the neighbour whose last binary digit is even, which is the one above. So in DKE Python, as in Python, three tenths make 0.30000000000000004. The number you meant is not the number you stored. A program that needs tenths to stay tenths has to count in something it can hold exactly: whole cents, or whole tenths of a second.

The Patriot's clock did count in whole tenths. Its trouble began when it converted them.

The last digit

Most of the time the last digit is noise, and nobody reads it. It matters when it is carried: added to itself a million times, compared with a threshold, or used to decide which of two answers a system believes. Then an error made the same way every time stops being small, and an answer that depends on which machine computed it stops being an answer you can check. The stricter test does not make arithmetic exact, and it cannot make a stored tenth a tenth. What it removes is the error the calculation itself adds, so that the last digit is the right one for the numbers given.

Every statement here about DKE is one its published documents make, at revision 133: the Python Language Specification and the Python Reference. Where they change, they govern. Python's answers were reproduced on CPython 3.11.2 and CPython 3.12.16, both over the GNU C Library 2.36 (Debian 12), and every exact value was computed independently to at least 80 digits and rounded once.