The Subtraction That Loses Every Correct Digit at Once

Two values agree to nine decimal digits. Subtract them, and the nine digits

they agreed on are gone from the answer and the digits they disagreed on are

all that is left. That much is obvious and is not the problem. The problem is

that the digits they disagreed on were never right in the first place, because

each value was only ever stored to within half a gap, and the subtraction has

just promoted that half gap from a rounding nobody could see into the entire

result.

This is the one operation in the whole of machine arithmetic that turns a small

relative error into a large one. Every other operation is tame. Multiply two

values each carrying a relative error of about $10^{-16}$ and the product

carries about $2 \times 10^{-16}$. Divide, take a square root, take a cosine,

and the same holds. Subtract two values that are nearly equal and the relative

error can be anything at all, including one hundred per cent.

Write $u$ for the unit roundoff of the 64 bit form, which is $2^{-53}$, or

0.00000000000000011102230246251565. That is the largest relative error a single

correctly rounded operation can leave behind, because the gap between

neighbouring values at magnitude $m$ is at most $m \cdot 2^{-52}$ and rounding

to nearest lands within half of it.

Let $a$ and $b$ be the exact values we meant, and $\hat{a}$ and $\hat{b}$ the

values the machine actually holds, so that $\hat{a} = a(1 + \delta_a)$ and

$\hat{b} = b(1 + \delta_b)$ with $|\delta_a| \le u$ and $|\delta_b| \le u$.

Here is the step that almost every explanation of this subject skips, and it is

the step that makes the rest make sense. If $\hat{a}$ and $\hat{b}$ lie within

a factor of two of each other, their difference is computed exactly. No

rounding happens. This is Sterbenz's lemma, and the reason is that the

difference of two such values needs no more significant bits than the smaller

of them already has, so the pattern for it already exists and the subtraction

simply writes it down.

So the computed difference is not an approximation of $\hat{a} - \hat{b}$. It

is $\hat{a} - \hat{b}$, which is

$$\hat{a} - \hat{b} = (a - b) + (a \delta_a - b \delta_b).$$

The absolute error in the difference is $|a \delta_a - b \delta_b|$, and it is

exactly the absolute error the operands brought with them. Nothing was added

and nothing was taken away. What changed is what that absolute error is being

compared against, because $a - b$ is tiny and $a$ and $b$ are not.