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.
\nFamously, with IEEE\ndouble-precision floats, 0.1 + 0.2 is not\n0.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.
\nIf 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:
\nWhere does this pattern come from?
\nA 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.
\nThe 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}})$.
\nThe effective value of the last fraction bit is thus $2^{E-52}$, a quantity called the\nulp, for “unit in the last place”, of $x$. The two floating-point numbers closest\nto $x$ differ from it by exactly one\nulp, except for an edge case on one side when $F = 0$ and $x$ is exactly a power of two. Example:
\n"hidden bit"52 bitsfloat(0.1) = 0.00011001100110011001100110011001100110011001100110011010₂ ulp(float(0.1)) = 0.00000000000000000000000000000000000000000000000000000001₂
\nThis 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.
\nAs an example (following e.g. qntm) let’s step through what\nhappens when you type 0.1 + 0.2 into the Python REPL.<sup>1</sup>
First, the Python expression 0.1 evaluates to some\nfloating-point number: specifically, as required by the IEEE standard,\nthe nearest floating-point number to 0.1. We’ll write that\nexact number as $\\text{float}(0.1)$.<sup>2</sup>
Similarly, the Python expression 0.2 evaluates to some\nother floating-point number, $\\text{float}(0.2)$.
Now, we can imagine the + being evaluated in two steps.\nFirst, Python computes the exact sum $\\text{float}(0.1) + \\text{float}(0.2)$.\nSecond, it rounds this to the nearest floating-point number. The\nresulting value is (\\text{float}(\\text{float}(0.1) +\n\\text{float}(0.2))). (This isn’t how it literally works — the\nexact sum from the first step isn’t ever materialized anywhere — but\nit’s mathematically accurate.)
Actually, there’s a subtlety here I did not notice until working this\nout in excruciating detail: $\\text{float}(0.1) + \\text{float}(0.2)$ is\nequally close to its two nearest floating-point numbers! When\nthis happens, Python rounds to the floating-point number with an even\nsignificand, in this case up.<sup>3</sup>
\nFinally, to print this value, Python has to convert it to a decimal.\nThis conversion is surprisingly subtle and not exactly specified by the\nIEEE standard! Even though (\\text{float}(\\text{float}(0.1) +\n\\text{float}(0.2))) isn’t exactly 0.3, it is pretty close, so it\nwouldn’t be unreasonable to display it as “0.3”. Another option would be\nto display it exactly, as\n“0.3000000000000000444089209850062616169452667236328125”. One might also\nimagine displaying it after various amounts of rounding:\n0.300000000000000044, or 0.30000000000000004441, or so on. One might\neven consider, say, 0.30000000000000005, because it is still true that\n(\\text{float}(\\text{float}(0.1) +\n\\text{float}(0.2)) = \\text{float}(0.30000000000000005)); that\nis, the Python expression 0.1 + 0.2 == 0.30000000000000005\nis true. The subtlety of this conversion is evidenced by the fact that,\nuntil a specific patch in\nPython 3.1, if you typed 1.1 into the REPL, Python\nwould print your input back to you as\n1.1000000000000001.
The standard description of how this decimal conversion should work\nwas formalized by Steele and White,\n1990<sup>4</sup>, who lay down three criteria<sup>5</sup>:
\nFollowing this algorithm, we can understand why\n0.1 + 0.2 == 0.30000000000000004.
Let’s do the general case: suppose you’re trying to add the exact\npositive decimal quantities $a$ and\n$b$, whose exact sum is $c$. Well, not fully general. We will\nassume that these values are “reasonable dollar amounts”, positive and\nless than $70 trillion; I think that should be enough to cover the\nreceipts I have to file. A bit above that (2<sup>46</sup> =\n70,368,744,177,664) floating-point numbers become sparser than multiples\nof cents, which is no good.
\nThe question is, how does $\\text{float}(\\text{float}(a) + \\text{float}(b))$ compare to $\\text{float}(c)$?
\nLet their difference be \n$$\n\\Delta := \\text{float}(\\text{float}(a) + \\text{float}(b)) - \\text{float}(c).\n$$\n We can reason about it by introducing the error function $\\text{error}(x) := \\text{float}(x) - x$. Then, we can rewrite $\\Delta$ as \n$$\n\\begin{aligned}\\Delta ={} &\\text{error}(a) + \\text{error}(b) \\&+ \\text{error}(\\text{float}(a) + \\text{float}(b)) - \\text{error}(c).\\qquad(*)\\end{aligned}\n$$
\nWe can bound each error term. Recall that the ulp (“unit in\nthe last place”) of a floating-point number is the value of the last\nbit. By mild abuse of notation, we will allow ourselves to write $\\text{ulp}(x)$ even when $x$ can’t be exactly represented as a\nfloating-point number, and understand that this means $\\text{ulp}(\\text{float}(x))$. So $\\text{float}(x) \\pm \\text{ulp}(x)$ are\nalso floating-point numbers<sup>6</sup>, which must be no closer\nto $x$ than $\\text{float}(x)$ itself (otherwise, $\\text{float}(x)$ would have evaluated to\nthe closer value); which means that, for all (reasonable) $x$, we have \n$$\n|\\text{error}(x)| \\leq\n\\frac{\\text{ulp}(\\text{float}(x))}{2}.\n$$\n Furthermore, equality\nonly holds when $x$ is exactly\nhalfway between the two closest floating-point numbers, which can’t hold\nif $x$ is a reasonable amount of\nmoney.<sup>7</sup> We can apply this bound term-by-term\nto $(*)$ to conclude that $|\\Delta| < 2\\text{ulp}(c)$.\nFurthermore, because $\\Delta$ is the\ndifference between two floating-point numbers near $c$, it’s a multiple of $\\text{ulp}(c)$.<sup>8</sup>\nFrom this we conclude that (\\Delta \\in\n{-\\text{ulp}(c), 0, +\\text{ulp}(c)}) — that is, the result can\nbe at most 1 ulp off from the answer.
\nHowever, 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:
\n$$\n\\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}\n$$\n Therefore, \n$$\n\\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}\n$$\n Subtracting $a + b = c$ from this, we get \n$$\n\\text{error}(a) + \\text{error}(b) \\leq \\text{error}(c) + \\text{ulp}(c).\n$$
\nThe same bound applies from the other side. As a result, if we let \n$$\n\\delta := \\text{error}(a) + \\text{error}(b) - \\text{error}(c),\n$$\n we have \n$$\n-1 \\leq \\frac{\\delta}{\\text{ulp}(c)} \\leq 1.\n$$\n 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)}$.
\nIf we make a heatmap of $\\delta/\\text{ulp}(c)$, we see what might be described as a more continuous version of Figure 1:
\nWe 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:
\nAnd, 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).
\nThe only remaining question is, why does $\\text{error}(x)$ look like that?
\nFirst 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 \n$$\n\\text{error}(x + 0.01) \\equiv \\text{error}(x) - 0.01 \\bmod \\text{ulp}(x);\n$$\n 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.
\nLet’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 \n$$\n\\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}\n$$\n 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]$.
\nTo understand the pattern even better, we can work out that \n$$\n\\frac{0.01 \\bmod 2^{-n}}{2^{-n}} = 0.01 \\times 2^n \\bmod 1 = \\frac{2^n \\bmod 100}{100}.\n$$\n 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 \n$$\n\\frac{0.01 \\bmod \\text{ulp}(2)}{\\text{ulp}(2)} \\equiv 0.48\\bmod 1\n$$\n 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.
\n(This is probably more practical than the main post)
\nHow 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$.
\nThe conceptually simplest way is to use some kind of exact rational\narithmetic, like Python’s fractions.\nFor my initial explorations, I used my own Noulith (after\nhaphazardly bolting on a bunch of features and bugfixes to its\nrational 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”.
\n2Sum\n(Møller, 1965): From $a$ and $b$, compute $s$ and $t$ such that (a\n+_\\text{float} b = s) and (a + b = s
\ndef two_sum(a: float, b: float) -> tuple[float, float]:\n s = a + b\n bb = s - a\n return s, (a - (s - bb)) + (b - bb)Veltkamp splitting<sup>9</sup>:\nFrom $a$, compute $h$ and $\\ell$ such that $a = h + \\ell$ exactly and both $h$ and $\\ell$ have at most 26 significant bits\n(after the hidden bit). This is useful because multiplying two such\nfloating-point numbers in floating-point is exact. (You can reallocate\nthe number of significant bits between $h$ and $\\ell$ by changing the magic constant.)
\ndef veltkamp_split(a: float) -> tuple[float, float]:\n c = ((1 << 27) | 1) * a\n hi = c - (c - a)\n return hi, a - hiDekker product<sup>10</sup>: From $a$ and $b$, compute $p$ and $r$ such that (a\n\\times_\\text{float} b = p) and (a\n\\times b = p + r) exactly.
\ndef dekker_product(a: float, b: float) -> tuple[float, float]:\n p = a * b\n ah, al = veltkamp_split(a)\n bh, bl = veltkamp_split(b)\n return p, ((ah * bh - p) + ah * bl + al * bh) + al * blUsing these techniques, we can compute a good-enough approximation to $\\text{error}(n / 100)$ as follows:
\ndef err_over_100(n: float) -> float:\n d = n / 100\n p, e = dekker_product(100, d)\n return ((p - n) + e) / 100