Exercises — Floating point gotchas — catastrophic cancellation, associativity failure

1 → (A). Both square roots are $\approx 10^6$ and nearly equal; subtracting near-twins reveals hidden relative error. Amplifier $\frac{|a|+|b|}{|a-b|}$ is huge. 2 → (B) absorption. Each term $10^{-9}$ is astronomically smaller than the running sum $\sim 10^9$. The gap between representable doubles near $10^9$ is about $10^9{\epsilon }_{\text{mach}}\approx 2\times 10^{-7}$, and $10^{-9}$ is far below that gap → each add contributes nothing. See Catastrophic cancellation's cousin, absorption. 3 → (C). Multiplication does not amplify relative error; the product is stored with relative error $\le u$. Safe. 4 → (A). $e^{10^{-10}}\approx 1+10^{-10}$; storing it rounds away almost all of the interesting part, then subtracting $1$ leaves a mangled result. Fix: use math.expm1.

1 → True. $\text{fl}(a+b)$ and $\text{fl}(b+a)$ round the same real number, so they agree exactly. 2 → False. Each + rounds at a place set by the partial sum's magnitude; regrouping rounds differently. 3 → False. By the Sterbenz lemma (defined at the top of this page), when $a,b$ are within a factor of $2$ the subtraction is exact. It merely exposes error already baked into the stored operands.

WHAT: plug numbers into the derived amplifier. $|a|+|b|=1.9999999$, and $|a-b|=10^{-7}$. $\frac{|a|+|b|}{|a-b|}=\frac{1.9999999}{10^{-7}}=1.9999999\times 10^7\approx 2.0\times 10^7.$ Worst-case relative error $=u\times 2.0\times 10^7\approx 1.11\times 10^{-16}\times 2.0\times 10^7\approx 2.2\times 10^{-9}$. WHY it matters: we started with $\sim 10^{-16}$ relative error and blew it up to $\sim 10^{-9}$ — we lost about 7 decimal digits, exactly ${\log }_{10}$ of the amplifier.

First expression: $1.0+10^{16}$. The true sum is $10^{16}+1$. The two candidate representable results are $10^{16}$ (below) and $10^{16}+2$ (above, one full gap up). The true value sits at $10^{16}+1$, which is exactly half a gap between them. IEEE-754 uses round-to-nearest, ties-to-even: it picks the neighbour whose last mantissa bit is $0$ (even). Here $10^{16}$ has last bit $0$ while $10^{16}+2$ has last bit... also representable, but $10^{16}$ is the even-mantissa neighbour, so the tie rounds down to $10^{16}$. Then $-10^{16}=0.0$. Prints 0.0. Second expression: $10^{16}-10^{16}=0.0$ exactly, then $+1.0=1.0$. Prints 1.0. Same three numbers, brackets only — answers differ by 100%. This is associativity failure via absorption. (The tie-to-even detail is why the $1$ vanishes rather than nudging up; even off the exact tie, any addend below $1$ would round down regardless.)

WHAT/WHY: $b<0$, so $-b>0$; the + branch adds same-sign numbers → safe root first. $\sqrt{b^2-4ac}=\sqrt{10^{16}-4}\approx 99999999.9999999998$, i.e. just under $10^8$. $x_{+}=\frac{-b+\sqrt{b^2-4ac}}{2a}=\frac{10^8+99999999.9999999998}{2}\approx \mathrm{99999999.99999999.}$ Now the dangerous root without subtraction: $x_{-}=\frac{c}{a{x}_{+}}=\frac{1}{99999999.99999999}\approx 1.0000000000000001\times 10^{-8}.$ Check: $x_{+}{x}_{-}=c/a=1$ ✓ and $x_{+}+x_{-}=-b/a=10^8$ ✓. The naive $x_{-}=\frac{-b-\sqrt{\dots }}{2a}=\frac{10^8-99999999.9999999998}{2}$ subtracts two $\sim 10^8$ numbers → catastrophic. See Quadratic formula numerical issues.

Naive: $\sqrt{10^{12}+1}\approx 1000000.0000005$ and $\sqrt{10^{12}}=1000000.0$. Their difference is $\sim 5\times 10^{-7}$, but each root has $\sim 10^{-16}$ relative error $\Rightarrow$ absolute error $\sim 10^{-10}$ in each; subtracting near-twins with amplifier $\frac{2\times 10^6}{5\times 10^{-7}}\approx 4\times 10^{12}$ wrecks it. Rewrite (WHY): multiply by the conjugate — this replaces subtraction of near-equals with an addition in the denominator: $g(x)=(\sqrt{x+1}-\sqrt{x})\cdot \frac{\sqrt{x+1}+\sqrt{x}}{\sqrt{x+1}+\sqrt{x}}=\frac{1}{\sqrt{x+1}+\sqrt{x}}.$ Stable value: $\frac{1}{\sqrt{10^{12}+1}+\sqrt{10^{12}}}\approx \frac{1}{2000000.0000005}\approx 5.0\times 10^{-7}$. Look at the figure: the red "near-twin" bars subtract to a sliver; the blue rewrite computes that sliver directly with no cancellation.

Why $0$: $\cos (10^{-8})=1-\frac{1}{2}(10^{-8}{)}^2+\cdots =1-5\times 10^{-17}$. But near $1.0$ the smallest step is ${\epsilon }_{\text{mach}}=2.22\times 10^{-16}$, and $5\times 10^{-17}$ is below half a step $\Rightarrow$ $\cos x$ rounds to exactly $1.0$. Then $1-1=0$, and $0/x^2=0$. Total loss. Stable: with $x/2=5\times 10^{-9}$, $\sin (x/2)/(x/2)\to 1$ (computed accurately for tiny arguments), so $\frac{1}{2}{\left(\frac{\sin (x/2)}{x/2}\right)}^2\approx \frac{1}{2}(1{)}^2=\mathrm{0.5.}$ Relative error $\sim u$ — perfect.

Rule: digits lost $\approx {\log }_{10}(\text{amplifier})$. ${\log }_{10}(10^9)=9$ digits lost. Survivors $\approx 16-9=7$ significant decimal digits. This is the practical meaning of the amplifier: its base-10 log is the digit death toll. See Round-off error propagation.

Recall the loop: y = x - c; t = s + y; c = (t - s) - y; s = t. Facts we use: the gap near $10^{16}$ is $2.0$, so integers there are only representable if even; $10^{16}$ and $10^{16}+2$ are representable, $10^{16}+1$ is not.

Naive left-to-right: $((1+10^{16})+1)-10^{16}$. $\text{fl}(1+10^{16})=10^{16}$ (the $1$ absorbed), $\text{fl}(10^{16}+1)=10^{16}$ (absorbed again), then $-10^{16}=0$. Naive result $=0.0$ — wrong by $2$.

Kahan trace (start $s=0.0,c=0.0$). Each line: WHAT the operation is, then the rounded value:

• $x=1$: y = 1 - 0 = 1 (exact). t = fl(0 + 1) = 1 (exact). c = fl((1 - 0) - 1) = 0. s = 1. → $s=1,c=0$.
• $x=10^{16}$: y = fl(1e16 - 0) = 1e16 (exact). t = fl(1 + 1e16) = 1e16 — the $1$ is below half the gap, absorbed. Now recover: t - s = fl(1e16 - 1) = 1e16 (again absorbed, exact difference $1e16-1$ rounds up to $1e16$), then c = fl((t - s) - y) = fl(1e16 - 1e16) = 0. Hmm — so on this simple integer data Kahan stores c = 0. → $s=10^{16},c=0$. (WHY the lost $1$ is not yet captured: both t-s and the recovery landed on exact grid points, so no residual surfaced this step — the loss is carried implicitly in the fact that s is $1$ short.)
• $x=1$: y = fl(1 - 0) = 1 (exact). t = fl(1e16 + 1) = 1e16 — absorbed again. t - s = fl(1e16 - 1e16) = 0, so c = fl(0 - 1) = -1. Now the dropped $1$ is captured: $c=-1$. s = 1e16. → $s=10^{16},c=-1$.
• $x=-10^{16}$: y = fl(-1e16 - (-1)) = fl(-1e16 + 1) = -1e16 (the $+1$ absorbed, but it is held in y's intent). t = fl(1e16 + (-1e16)) = 0 (exact cancellation). t - s = fl(0 - 1e16) = -1e16, c = fl((-1e16) - (-1e16)) = 0. s = 0. → $s=0,c=0$.

Final s = 0.0 with this exact IEEE trace. The single interleaving of $\pm 10^{16}$ around lone $1$'s defeats even Kahan here, because the compensation $c=-1$ is itself absorbed when re-added to $-10^{16}$. Lesson: Kahan is not magic against pure absorption of a subnormal-relative addend into a giant that is then cancelled — the cure is to reorder (cancel the giants first). This is exactly Exercise 5.3.

Contrast — an order Kahan does save. On $[10^{16},1,1,-10^{16}]$, naive gives $0$, but summing $[1,1,10^{16},-10^{16}]$ left-to-right gives $\text{fl}(\text{fl}(\text{fl}(1+1)+10^{16})-10^{16})=\text{fl}(\text{fl}(2+10^{16})-10^{16})$; since $10^{16}+2$ is representable, this is $\text{fl}((10^{16}+2)-10^{16})=2$. Exact. The moral holds: order beats absorption; Kahan helps when losses are residuals, not total absorptions of terms that later cancel.

Predicate: $|x-y|\le \max (\text{rtol}\cdot \max (|x|,|y|),\text{atol}).$

• rtol (relative tolerance) scales with the operands' magnitude — needed because float error is relative.
• atol (absolute tolerance) rescues the case where both are near $0$ (relative test degenerates there). The $0.1+0.2$ case: the computed sum is $0.30000000000000004$, differing from $0.3$ by $\approx 5.55\times 10^{-17}$. Threshold $=1\times 10^{-9}\times 0.3=3\times 10^{-10}$. Since $5.55\times 10^{-17}\le 3\times 10^{-10}$, they ARE close → predicate returns True, even though == returns False.

Step 1 (WHAT): expand the computed result. $\hat{r}=\hat{a}-\hat{b}=a(1+{\delta }_a)-b(1+{\delta }_b)=(a-b)+(a{\delta }_a-b{\delta }_b).$ Step 2: subtract the true $r=a-b$: $\hat{r}-r=a{\delta }_a-b{\delta }_b.$ Step 3 (bound the numerator): by the triangle inequality and $|\delta |\le u$, $|a{\delta }_a-b{\delta }_b|\le |a||{\delta }_a|+|b||{\delta }_b|\le u(|a|+|b|).$ Step 4 (divide by $|r|=|a-b|$): $\left|\frac{\hat{r}-r}{r}\right|=\frac{|a{\delta }_a-b{\delta }_b|}{|a-b|}\le u\frac{|a|+|b|}{|a-b|}.■$ WHY the "exact subtraction" assumption is fair: the Sterbenz lemma guarantees no new rounding when $\frac{1}{2}a\le b\le 2a$; all the damage lives in ${\delta }_a,{\delta }_b$, which is precisely what the proof isolates.

Step 1 (WHAT the gap is): near $L$ the spacing between representable doubles is $\text{gap}(L)=2^{\lfloor {\log }_{2}L\rfloor }\cdot 2^{-52}=2^{\lfloor {\log }_{2}L\rfloor -52}.$ This is the weight of the least-significant mantissa bit at $L$'s exponent. Step 2 (round-to-nearest rule): the true value $L+s$ lies between the two neighbours $L$ and $L+\text{gap}(L)$. Round-to-nearest snaps it back to $L$ whenever it is at most half a gap above $L$: $s\le \frac{1}{2}\text{gap}(L)=2^{\lfloor {\log }_{2}L\rfloor -53}.$ At exactly half a gap, ties-to-even decides: if $L$'s last mantissa bit is $0$ (even) the tie also rounds to $L$, so the strict absorption region is $s\le \frac{1}{2}\text{gap}(L)$ for even-$L$ (and $s<\frac{1}{2}\text{gap}(L)$ otherwise). Step 3 (apply $L=2^{53}$): $\lfloor {\log }_{2}L\rfloor =53$, so $\text{gap}(2^{53})=2^{53-52}=2^1=2,\frac{1}{2}\text{gap}=2^{53-53}=2^0=1.$ Therefore any $s\le 1$ is fully absorbed: $\text{fl}(2^{53}+1)=2^{53}$. Indeed $2^{53}+1$ is the first positive integer a double cannot represent — a famous landmark: doubles represent every integer up to $2^{53}$, then start skipping.

Bad order $[10^{16},1,1,-10^{16}]$: absorb, absorb, then cancel → $0$. Good order $[10^{16},-10^{16},1,1]$: cancel the giants first ($10^{16}-10^{16}=0$ exactly), then $0+1+1=2$. Exact. Principle: kill large opposite-sign terms against each other before they can swallow the small survivors. When you can't reorder (streaming data), use Kahan for residual losses — but as Exercise 4.1 showed, even Kahan cannot rescue a total absorption that is later cancelled, so reordering is the primary weapon. This is the associativity-failure lesson turned into a strategy.

• Digits lost from an amplifier of $10^9$? ::: about $9$ (i.e. ${\log }_{10}$ of it), leaving $\sim 7$.
• Stable rewrite of $\sqrt{x+1}-\sqrt{x}$? ::: $\frac{1}{\sqrt{x+1}+\sqrt{x}}$.
• First integer a double cannot represent? ::: $2^{53}+1$.
• math.isclose predicate? ::: $|x-y|\le \max (\text{rtol}\cdot \max (|x|,|y|),\text{atol})$.
• Exact sum of $[10^{16},-10^{16},1,1]$ left-to-right? ::: $2.0$ (cancel giants first).
• Which quadratic root do you compute directly, sign-wise? ::: the one whose $\pm$ matches the sign of $-b$ (same-sign addition).