Adding Floating-Point Decimals for Fun and Profit

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.¹

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)$.²

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.³

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, 1990 ⁴, who lay down three criteria⁵:

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 (2⁴⁶ = 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 numbers⁶, 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.⁷ 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)$.⁸ 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) - 0.01 \bmod \text{ulp}(x);\] that is, $\text{error}(x)$ is an “arithmetic sequence with common difference −0.01” modulo $\text{ulp}(x)$. So, the wraparound behavior of this function leads to the periodic patterns in our previous figures.

Let’s focus on the lower-right quadrant of Figure 2, $[0.5, 1] \times [0.5, 1]$. In this region we can compute that $\text{ulp}(0.5) = 2^{-53}$ and $\text{ulp}(1) = 2^{-52}$, and then that \[\begin{aligned}\frac{0.01 \bmod \text{ulp}(0.5)}{\text{ulp}(0.5)} &\equiv 0.92\equiv -0.08 \bmod 1 \\ \frac{0.01 \bmod \text{ulp}(1)}{\text{ulp}(1)} &\equiv 0.96\equiv -0.04 \bmod 1,\end{aligned}\] which are both “close to 0”. Because $0.08 \approx 1/12$, $\text{error}(x)$ has 12 “steps” before wrapping around when $x \in [0.5, 1]$; and because $0.04 = 1/25$, $\text{error}(x)$ has 25 “steps” before wrapping around when $x \in [1, 2]$.

To understand the pattern even better, we can work out that \[\frac{0.01 \bmod 2^{-n}}{2^{-n}} = 0.01 \times 2^n \bmod 1 = \frac{2^n \bmod 100}{100}.\] It is actually a nice coincidence that the number of fraction bits in double-precision floating-point, 52, is such that $2^{52}$ is “close to 0” mod 100; that’s the reason the error function doesn’t wrap around so much, so we have smooth regions. If we expand our diagrams to $a, b \in [0, 2]$ such that $c$ can reach $[2, 4]$, we see messier checkerboards and diagonal lines, because \[\frac{0.01 \bmod \text{ulp}(2)}{\text{ulp}(2)} \equiv 0.48\bmod 1\] and $\text{error}(x)$ wraps around roughly every other step in $[2, 4]$, which then interferes with the parity of $c$’s significand in a more complicated way.

Figure 6: How floating-point affects summing two multiples of 0.01, up to 2.00

Figure 7: Combined floating-point error from summing two multiples of 0.01, up to 2.00

Figure 8: Multiples of 0.01 with odd floating-point significand, up to 4.00

Figure 9: The error function at multiples of 0.01, up to 4.00

Appendix: Error-free transformations in floating-point

(This is probably more practical than the main post)

How do you actually calculate a function like $\text{error}(x) = \text{float}(x) - x$ on a computer, for example, to generate the figures in this post? Obviously you can’t directly compute it in the same floating-point format you’re studying. In that format, $\text{float}$ is the identity function; the error has already been incurred by the time you try to express $x$.

The conceptually simplest way is to use some kind of exact rational arithmetic, like Python’s fractions. For my initial explorations, I used my own Noulith (after haphazardly bolting on a bunch of features and bugfixes to its rational type…).

However, it turns out there are a bunch of indirect ways to work with errors like this without leaving the floating-point format. I believe these techniques are called “error-free transformations”.

2Sum (Møller, 1965): From $a$ and $b$, compute $s$ and $t$ such that $a +_\text{float} b = s$ and $a + b = s + t$ exactly.

def two_sum(a: float, b: float) -> tuple[float, float]:
    s = a + b
    bb = s - a
    return s, (a - (s - bb)) + (b - bb)

Veltkamp splitting⁹: From $a$, compute $h$ and $\ell$ such that $a = h + \ell$ exactly and both $h$ and $\ell$ have at most 26 significant bits (after the hidden bit). This is useful because multiplying two such floating-point numbers in floating-point is exact. (You can reallocate the number of significant bits between $h$ and $\ell$ by changing the magic constant.)

def veltkamp_split(a: float) -> tuple[float, float]:
    c = ((1 << 27) | 1) * a
    hi = c - (c - a)
    return hi, a - hi

Dekker product¹⁰: From $a$ and $b$, compute $p$ and $r$ such that $a \times_\text{float} b = p$ and $a \times b = p + r$ exactly.

def dekker_product(a: float, b: float) -> tuple[float, float]:
    p = a * b
    ah, al = veltkamp_split(a)
    bh, bl = veltkamp_split(b)
    return p, ((ah * bh - p) + ah * bl + al * bh) + al * bl

Using these techniques, we can compute a good-enough approximation to $\text{error}(n / 100)$ as follows:

def err_over_100(n: float) -> float:
    d = n / 100
    p, e = dekker_product(100, d)
    return ((p - n) + e) / 100

I’ll be saying “Python” to be concrete and concise, but most of this process is the same in any programming language or computing environment with IEEE double-precision floats (which is a lot of them), since they adhere to the same standard and probably rely on the same CPU instructions. In particular, JavaScript, which is computing the figures in this post live in your browser. (This is a somewhat irresponsible way to generate them when they’re all theoretically static, I’m sorry)↩︎
If we wanted, we could write $\text{float}(0.1)$ out as an exact decimal: $\text{float}(0.1) = 0.1000000000000000055511151231257827021181583404541015625$. This is neat but a bit unwieldy.↩︎

This rounding strategy is customizable, if you have fairly low-level access to your FPU. Here’s a C program that demonstrates it with fesetround. (I use an obscenely high precision in the printf because changing the rounding mode also affects how the float gets printed as a decimal, but that’s a red herring.)

#pragma STDC FENV_ACCESS on

#include <stdio.h>
#include <fenv.h>

int main() {
    double a = 0.1, b = 0.2;
    printf("0.3         :  %.54g\n", 0.3);
    printf("FE_TONEAREST:  %.54g\n", a + b);
    fesetround(FE_UPWARD);
    printf("FE_UPWARD:     %.54g\n", a + b);
    fesetround(FE_DOWNWARD);
    printf("FE_DOWNWARD:   %.54g\n", a + b);
    fesetround(FE_TOWARDZERO);
    printf("FE_TOWARDZERO: %.54g\n", a + b);
}

I don’t think Python directly exposes this as a setting, but you can muck around with ctypes to get the same effect:

import ctypes
import ctypes.util

libm = ctypes.CDLL(ctypes.util.find_library('m'))

x, y = 0.1, 0.2
print(x + y)
libm.fesetround(0x400) # FE_DOWNWARD, from fenv.h
print(x + y)

These should bottom out into instructions like x86’s FLDCW and LDMXCSR.↩︎

though they note that it was folkloric knowledge before they published↩︎
They have a fourth criterion about the algorithm for generating the digits, namely that it should do so from left to right without backtracking. We won’t concern ourselves with that here.↩︎
In the edge case where $x$ is a power of two, this is still true, although $\text{float}(x) - \text{ulp}(\text{float}(x))$ won’t be one of the closest floating-point numbers to $x$.↩︎
The smallest multiple of 0.01 that’s exactly halfway between two representable floating-point numbers is around 2 quadrillion: 2,251,799,813,685,248.25 = 2⁵¹ + ¼.↩︎
I’m ignoring the case where $c$ is a power of two again. It’s annoying and left as an exercise to the reader.↩︎
G. W. Veltkamp, “ALGOL procedures voor het berekenen van een inwendig product in dubbele precisie,” RC-Informatie, Technishe Hogeschool Eindhoven, Tech. Rep. 22, 1968.↩︎
T. J. Dekker, “A floating-point technique for extending the available precision,” Numerische Mathematik, vol. 18, no. 3, pp. 224–242, 1971.↩︎