Adding Floating-Point Decimals for Fun and Profit

2026-08-30 (2527 words) filed under Math, CS Many people know that you shouldn’t do decimal calculations, such as those involving U.S. dollars and cents, with the floating-point numbers in most programming languages. This is because decimal numbers can’t be expressed exactly as such floating-point numbers, so you will encounter rounding errors. Famously, with IEEE double-precision floats, 0.1 + 0.2 is not 0.3, but rather, 0.30000000000000004. On the other hand, I use a Python REPL to add up decimal numbers for receipts all the time. Of course, my stakes are lower; I’m summing things in the $1–$100 range and know to manually round the microscopic errors off before copying the sum somewhere. But actually it’s quite often that those errors don’t appear at all. If we sum every pair of multiples of 0.01 up to 1.00 and look at whether the printed result is too large (blue), too small (red), or correct, we get a cool pattern: Figure 1: How floating-point affects summing two multiples of 0.01, up to 1.00 Where does this pattern come from? Floats, briefly A double-precision floating-point number \(x\) consists of 1 sign bit, 11 exponent bits, and 52 fraction bits (in that order), for a total of 64 bits. The sign bit provides a sign, \(+\) or \(-\). The exponent bits represent an integer \(E\) between −1022 and 1023, inclusive. The fraction bits represent a nonnegative integer \(F\) less than \(2^{52}\), which corresponds to the significand \(1 + F/{2^{52}} \in [1, 2)\); that ever-present \(1\) is called the “hidden bit”. The floating-point number’s value is \(x = \pm 2^E(1 + F/{2^{52}})\). The effective value of the last fraction bit is thus \(2^{E-52}\), a quantity called the ulp, for “unit in the last place”, of \(x\). The two floating-point numbers closest to \(x\) differ from it by exactly one ulp, except for an edge case on one side when \(F = 0\) and \(x\) is exactly a power of two. Example: "hidden bit"52 bits float(0.1)  = 0.00011001100110011001100110011001100110011001100110011010₂ ulp(float(0.1)) = 0.00000000000000000000000000000000000000000000000000000001₂ This description is good enough for our purposes but ignores many other cases: 0, subnormal numbers, infinities, and NaNs; they use the two values of exponent bits I didn’t describe. I also won’t consider other precisions of floating-point numbers, for simplicity. Anatomy of an addition As an example (following e.g. qntm) let’s step through what happens when you type 0.1 + 0.2 into the Python REPL.1 First, the Python expression 0.1 evaluates to some floating-point number: specifically, as required by the IEEE standard, the nearest floating-point number to 0.1. We’ll write that exact number as \(\text{float}(0.1)\).2 Similarly, the Python expression 0.2 evaluates to some other floating-point number, \(\text{float}(0.2)\). Now, we can imagine the + being evaluated in two steps. First, Python computes the exact sum \(\text{float}(0.1) + \text{float}(0.2)\). Second, it rounds this to the nearest floating-point number. The resulting value is \(\text{float}(\text{float}(0.1) + \text{float}(0.2))\). (This isn’t how it literally works — the exact sum from the first step isn’t ever materialized anywhere — but it’s mathematically accurate.) Actually, there’s a subtlety here I did not notice until working this out in excruciating detail: \(\text{float}(0.1) + \text{float}(0.2)\) is equally close to its two nearest floating-point numbers! When this happens, Python rounds to the floating-point number with an even significand, in this case up.3 Finally, to print this value, Python has to convert it to a decimal. This conversion is surprisingly subtle and not exactly specified by the IEEE standard! Even though \(\text{float}(\text{float}(0.1) + \text{float}(0.2))\) isn’t exactly 0.3, it is pretty close, so it wouldn’t be unreasonable to display it as “0.3”. Another option would be to display it exactly, as “0.3000000000000000444089209850062616169452667236328125”. One might also imagine displaying it after various amounts of rounding: 0.300000000000000044, or 0.30000000000000004441, or so on. One might even consider, say, 0.30000000000000005, because it is still true that \(\text{float}(\text{float}(0.1) + \text{float}(0.2)) = \text{float}(0.30000000000000005)\); that is, the Python expression 0.1 + 0.2 == 0.30000000000000005 is true. The subtlety of this conversion is evidenced by the fact that, until a specific patch in Python 3.1, if you typed 1.1 into the REPL, Python would print your input back to you as 1.1000000000000001. The standard description of how this decimal conversion should work was formalized by Steele and White, 19904, who lay down three criteria5: The decimal representation should round-trip: if you type it in again, you should get the same floating-point number. This rules out the output “0.3”. Subject to criterion 1, the decimal representation should be the shortest possible. This rules out outputs like “0.30000000000000004441”. Subject to criteria 1 and 2, the decimal representation should be as close as possible to the floating-point number. This rules out outputs like “0.30000000000000005”. Following this algorithm, we can understand why 0.1 + 0.2 == 0.30000000000000004. Generalizing Let’s do the general case: suppose you’re trying to add the exact positive decimal quantities \(a\) and \(b\), whose exact sum is \(c\). Well, not fully general. We will assume that these values are “reasonable dollar amounts”, positive and less than $70 trillion; I think that should be enough to cover the receipts I have to file. A bit above that (246 = 70,368,744,177,664) floating-point numbers become sparser than multiples of cents, which is no good. The question is, how does \(\text{float}(\text{float}(a) + \text{float}(b))\) compare to \(\text{float}(c)\)? Let their difference be \[\Delta := \text{float}(\text{float}(a) + \text{float}(b)) - \text{float}(c).\] We can reason about it by introducing the error function \(\text{error}(x) := \text{float}(x) - x\). Then, we can rewrite \(\Delta\) as \[\begin{aligned}\Delta ={} &\text{error}(a) + \text{error}(b) \\&+ \text{error}(\text{float}(a) + \text{float}(b)) - \text{error}(c).\qquad(*)\end{aligned}\] We can bound each error term. Recall that the ulp (“unit in the last place”) of a floating-point number is the value of the last bit. By mild abuse of notation, we will allow ourselves to write \(\text{ulp}(x)\) even when \(x\) can’t be exactly represented as a floating-point number, and understand that this means \(\text{ulp}(\text{float}(x))\). So \(\text{float}(x) \pm \text{ulp}(x)\) are also floating-point numbers6, which must be no closer to \(x\) than \(\text{float}(x)\) itself (otherwise, \(\text{float}(x)\) would have evaluated to the closer value); which means that, for all (reasonable) \(x\), we have \[|\text{error}(x)| \leq \frac{\text{ulp}(\text{float}(x))}{2}.\] Furthermore, equality only holds when \(x\) is exactly halfway between the two closest floating-point numbers, which can’t hold if \(x\) is a reasonable amount of money.7 We can apply this bound term-by-term to \((*)\) to conclude that \(|\Delta| < 2\text{ulp}(c)\). Furthermore, because \(\Delta\) is the difference between two floating-point numbers near \(c\), it’s a multiple of \(\text{ulp}(c)\).8 From this we conclude that \(\Delta \in \{-\text{ulp}(c), 0, +\text{ulp}(c)\}\) — that is, the result can be at most 1 ulp off from the answer. However, here’s a derivation that produces tighter intermediate bounds on \(|\Delta|\): Assume without loss of generality that \(a \leq b\). Then, \(\text{float}(c) + \text{ulp}(c) - \text{float}(b)\) is a representable floating-point number because the result’s ulp is ≤ that of both \(b\) and \(c\). Therefore, at least that is an available approximation of \(a\). And it’s an overestimate: \[\begin{aligned}a &= c - b \\ &\leq \text{float}(c) + \frac{\text{ulp}(c)}{2} - \text{float}(b) + \frac{\text{ulp}(c)}{2} \\ &= \text{float}(c) - \text{float}(b) + \text{ulp}(c).\end{aligned}\] Therefore, \[\begin{aligned}\text{float}(a) &\leq \text{float}(c) - \text{float}(b) + \text{ulp}(c)\\ \text{float}(a) + \text{float}(b) &\leq \text{float}(c) + \text{ulp}(c).\end{aligned}\] Subtracting \(a + b = c\) from this, we get \[\text{error}(a) + \text{error}(b) \leq \text{error}(c) + \text{ulp}(c).\] The same bound applies from the other side. As a result, if we let \[\delta := \text{error}(a) + \text{error}(b) - \text{error}(c),\] we have \[-1 \leq \frac{\delta}{\text{ulp}(c)} \leq 1.\] As before, we know \(|\delta - \Delta| \leq \text{ulp}(c)/2\), and again since \(\Delta\) is a multiple of \(\text{ulp}(c)\) we see that \(\Delta \in \{-\text{ulp}(c), 0, +\text{ulp}(c)\}\). If we make a heatmap of \(\delta/\text{ulp}(c)\), we see what might be described as a more continuous version of Figure 1: Figure 2: Combined floating-point error from summing two multiples of 0.01, up to 1.00 We can now understand Figure 1 as a “rounded” version of Figure 2, with a checkerboard pattern arising in regions where \(\delta\) is exactly \(\pm\text{ulp}(c)/2\) due to rounding to floats with even significand: Figure 3: Multiples of 0.01 with odd floating-point significand And, we can interpret Figure 2 as the result of “interference” between three copies of the \(\text{error}\) function: one horizontal, one vertical, one diagonal (albeit with a changing denominator). Figure 4: Decomposing Figure 2 into terms The only remaining question is, why does \(\text{error}(x)\) look like that? One-dimensional error Figure 5: The error function at multiples of 0.01, up to 2.00 First let’s observe that \(\text{error}(x) = 0\) whenever \(x\) is an exact power of 2. In between two such powers, let’s compare \(\text{error}(x)\) and \(\text{error}(x+0.01)\). We have \(\text{ulp}(x) = \text{ulp}(x+0.01)\), so \(\text{float}(x + 0.01) \equiv 0 \equiv \text{float}(x) \bmod \text{ulp}(x)\), so \[\text{error}(x + 0.01) \equiv \text{error}(x) -