Hi Friends,
Welcome to the 196th issue of the Polymathic Engineer newsletter. This time, we keep talking about floating-point numbers and see how the imprecision of floating-point arithmetic can have dangerous effects. On February 1991, during the first Gulf War, an American Patriot missile battery in Saudi Arabia, failed to intercept an incoming Iraqi Scud missile. The Scud struck an Army barracks and killed many soldiers. The US analyzed the failure and determined that the root cause was an imprecision in a numeric calculation.
The Patriot system had an internal clock implemented with a counter incremented every tenth of a second. To get the time in seconds, the software multiplied the counter by 0.1, stored in a 24-bit fixed-point register with 23 bits after the binary point. However, 0.1 has no finite binary expansion, so the register held a truncated value, short by about 0.000000095. That is a tenth of a microsecond, and nothing to worry about for one tick. But the battery had been running for about 100 hours, 3.6 million ticks, and by then the computed time was 0.34 seconds off. The tracking software used that clock to work out where the incoming missile would be, and it was estimated that the window it searched was off by 687 meters. The Scud was not where the Patriot looked.
Every operation involving floats computes the exact result and rounds it to the nearest value the format can represent. Here, we look at what happens when you compute with these numbers: how an error that is harmless in a single operation accumulates, and why comparing two floats can be disastrous.
The outline is as follows:
Adding floats
Multiplying floats
Comparing floats
Agents Are Joining the Engineering Team
The future of software engineering isn’t about replacing people with agents. It’s about giving agents well-defined jobs so teams can work better together. Join Honeycomb co-founders Charity and Christine, November 3-4, for Innovation Week: a two-day virtual event with live Q&A exploring what AI SRE looks like done right, from the agents on your team to visibility across your whole ecosystem.
Thanks to Honeycomb for making this issue free.
Adding floats
Let’s start with addition, using as an example a mantissa with 3 significant digits and a 2-digit exponent. Since the format supports signed values, subtraction works the same way as addition: X − Y is just X + (−Y).
Whenever you add or subtract two numbers, you must first adjust the two values so that they have the same exponent. For example, consider the addition of 1.23e1 and 4.56e0. To add these two numbers, you can increase the exponent of the second number by shifting its decimal point to the left, converting it to 0.456e1, then add them to get 1.686e1. The result does not fit into the three significant digits of our format, so you must round it to 1.69e1. As you can see, the lack of precision (the number of digits kept in a computation) affects the accuracy (the correctness of the computation).
In the previous example, we were able to round the result because we maintained four significant digits during the calculation. If our floating-point calculation were limited to three significant digits during computation, we would have had to truncate the last digit of the smaller number, obtaining 1.68e1, which is even less correct. To improve the accuracy of floating-point calculations, it is necessary to use extra digits during the calculation. These extra digits are known as guard digits (or guard bits in the case of a binary format). They greatly enhance accuracy during a long chain of computations.
The accuracy lost in a single operation is usually small. However, if you compute a value that is the result of a sequence of operations, the error accumulates and greatly affects the computation itself. Suppose we add 1.23e3 and 1.00e0. Adjusting the numbers so their exponents are the same before the addition produces 1.23e3 + 0.001e3 = 1.231e3, which rounds back to 1.23e3. The 1 disappeared entirely. This might seem reasonable; after all, if we only keep three digits, adding a small value shouldn’t affect the result. However, suppose we add 1.00e0 to 1.23e3 ten times. The first time we get 1.23e3, and we get the same result the second, third, and tenth time. Had we added 1.00e0 to itself ten times first, then added the result 1.00e1 to 1.23e3, we would get 1.24e3. The order of evaluation changes the result. The reason is the one we saw on the number line of the previous issue: representable values get sparser as they grow, and around 1,230 the gap between neighbors is wider than the 1 we are adding.
Floating-point addition is commutative, but it is not associative, and that has important implications for scientific programmers and compiler writers. A compiler can’t rewrite x = a + b + c and y = b + c + d to compute b + c once, since the two groupings may yield different values. As a result, compilers tend to be very conservative, avoiding any optimizations that could have even the slightest effect on functionality. And summing a list gives a different total depending on the order: in Python, sum([0.1] * 10) gives 0.9999999999999999, while math.fsum, which keeps track of the lost digits, gives 1.0.
Subtraction has a problem of its own. Consider the computation 1.23e0 − 1.22e0, which produces 0.01e0. Although this is mathematically equivalent to 1.00e−2, this form suggests that the last two digits are both exactly zero, when we only have a single significant digit left. Whenever you subtract two numbers that are close to each other, most of the digits cancel out, and the accuracy of the result may be far less than the precision available in the format. The digits that survive carry no information. This is known as cancellation, and it is one of the main ways a long computation goes wrong.
In general, your results will be better when adding or subtracting numbers whose relative magnitudes (that is, the sizes of the exponents) are similar. In a long chain of additions and subtractions, group the operations so that values close to one another in magnitude meet first, and keep an eye on any step that subtracts two nearly equal values.
Multiplying floats
Multiplication and division do not suffer from the alignment problem, since you do not have to adjust the exponents before the operation. All you need to do is add the exponents and multiply the mantissas (or subtract the exponents and divide the mantissas). By themselves, multiplication and division do not produce particularly poor results. However, they tend to make worse any error that already exists in a value. For example, if you multiply 1.23e0 by 2 when you should be multiplying 1.24e0 by 2, the result is twice as far from the truth as it was.
This is precisely what happened to the Patriot. The stored 0.1 was wrong by a tenth of a microsecond, and multiplying it by a counter of 3.6 million scaled that single error up to a third of a second. An offset like this would normally not matter, since tracking a target depends on the difference between two successive clock readings, and the same offset cancels out. The problem was that the Patriot software had been upgraded to use a more accurate function for reading the time, but not all of the function calls had been replaced by the new code. As a result, the tracking software used the accurate time for one reading and the inaccurate time for the other, and the difference between them was off by the whole third of a second.
Multiplication and division have another problem. When multiplying two very large or two very small numbers, it is quite possible for the result to overflow or underflow, that is, to go beyond the largest value the format can represent or below the smallest. The same happens when dividing a small number by a large number, or a large number by a small one. In double precision, (1e200 × 1e200) × 1e−200 evaluates to infinity, since the intermediate product overflows, while 1e200 × (1e200 × 1e−200) evaluates to 1e200. So multiplication is not associative either, and it does not distribute over addition: 1e200 × 1e200 − 1e200 × 1e200 is ∞ − ∞, which is NaN. The rule here is the same as for addition. When multiplying and dividing sets of numbers, try to pair numbers that have similar magnitudes, and order the operations so that intermediate results stay in range.
Converting between formats deserves a mention too, since it is a rounding step that is easy to forget. Going from an int to a float cannot overflow, but it may round: a float keeps 24 significant bits, so every integer up to 16,777,216 is exact, and 16,777,217 becomes 16,777,216. Going from an int or a float to a double preserves the exact value, since a double has both greater range and greater precision. Going from a double to a float, the value can overflow to +∞ or −∞, since the range is smaller, and otherwise it may be rounded. Going from a float or a double to an int, the value is rounded toward zero: 1.999 becomes 1 and −1.999 becomes −1. The value may also overflow, and the C standard does not specify a result for this case.
Comparing floats
Comparing floating-point numbers is very dangerous. Given the inaccuracies present in any computation, including converting the literal you typed into a floating-point value, you should never compare two floats to see if they are equal. Different computations that produce the same mathematical result may differ in their least significant bits, and the test for equality succeeds only if all the bits in the two operands are the same. In Python, 0.1 * 3 == 0.3 is False, since the product rounds to 0.30000000000000004 while the literal 0.3 rounds to 0.29999999999999999. The straight comparison fails when, algebraically, it should succeed.
The standard way to test for equality is to decide how much error you are willing to tolerate and then check whether one value is within that error of the other: abs(a − b) <= error. You must exercise care when choosing the value for error. It should be slightly greater than the largest amount of error that will creep into your computations, and that depends on the format you use and on the magnitudes of the values you compare. A fixed tolerance like 1e−9 is far too loose for values around 1e−12 and far too strict for values around 1e12. For this reason, the tolerance is usually expressed relative to the values themselves, as a fraction of the larger of the two: abs(a − b) <= rel × max(abs(a), abs(b)). Python’s math.isclose does this, with an optional absolute tolerance for comparisons near zero, where a relative one would shrink to nothing.
Checking two floats for equality is a famous problem, and almost every introductory programming text discusses it. Less well known is the fact that comparing for less than or greater than creates the same problems. Suppose that some chain of calculations produces the result 1.25, which is only accurate to plus or minus 0.05, so the real value could be anywhere between 1.20 and 1.30. A second chain of calculations produces 1.27, accurate to the full precision of the format. If you compare the two, you find that the first is less than the second. Given the inaccuracy of the first calculation, this might not be true: if its correct result was 1.28, the first is greater. The only reasonable test is to first check whether the two values are within the error tolerance of one another. If so, treat them as equal. Only if they are not, can you compare them to see which one is less than the other.
Wrapping Up
Floating-point arithmetic is an approximation of real arithmetic, which holds up well for a single operation. The trouble is a long chain of them, where the rounding meet: through alignment, when a small term is added to a big one; through cancellation, when two close values are subtracted; and through multiplication, when a tiny error in one operand is scaled by a big factor.
The Patriot’s third of a second was one truncation error of 0.1, less than a tenth of a microsecond, multiplied by 3.6 million ticks. Nothing in the software was broken, and every step did what the format said it would do.