The trap: subtracting two numbers that are almost equal
By now you know the rounded-operation model from earlier in this rung: every stored number and every basic operation carries a tiny relative error of size at most the unit roundoff u, about 10^(-16) in double precision. You might expect that as long as each step is accurate to 16 digits, the answer stays accurate to 16 digits. That hope breaks in exactly one common situation, and the break is dramatic. It has a name worth remembering: catastrophic cancellation.
Here is the whole idea in one tiny picture. Suppose two numbers agree in their first several digits — say a = 1.2345678 and b = 1.2345600. Each is fine on its own. But subtract them and you get a - b = 0.0000078. Seven matching leading digits vanished; the result keeps only the two trailing digits that differed. Those trailing digits are precisely the ones least protected by rounding — the noise you couldn't see before is suddenly the whole answer.
Notice the sleight of hand: the subtraction itself was computed almost perfectly. If a and b were the exact numbers you wanted, fl(a - b) is correct to within u, no problem. The damage was already done when a and b were stored or computed: each carried an error of relative size u, which in absolute terms is roughly 10^(-16) times their magnitude (~1.2). After cancellation the answer is only ~10^(-5), so that same absolute error is now a relative error of about 10^(-16)/10^(-5) = 10^(-11) — eleven of your sixteen digits, gone. Subtraction did not create the error; it merely promoted the existing rounding noise from invisible to dominant.
Conditioning vs. stability: whose fault is the error?
It helps to name two different things that can go wrong, because catastrophic cancellation is a stability issue, not a conditioning one. A problem's condition number measures how much the true answer reacts to tiny changes in the inputs — it belongs to the problem, no matter how you solve it. The numerical stability of an algorithm measures how much extra error the algorithm injects on top of that. The honest accounting rule from earlier: accuracy ≈ conditioning × stability. A backward-stable algorithm still gives a bad answer on an ill-conditioned problem — and nothing can save you there.
Catastrophic cancellation is the classic example of an unstable algorithm hurting a perfectly well-conditioned problem. Computing 0.0000078 from inputs near 1.23 is itself a benign, well-conditioned task — a small relative change in the true a and b would change the true difference by a similarly small relative amount. The instability lives entirely in the way we computed it: we routed the calculation through a subtraction that magnifies the inputs' rounding error. Same problem, different algorithm, hugely different forward error.
The quadratic formula, and where it falls apart
Now the famous victim. To solve a x^2 + b x + c = 0 you reach, automatically, for the formula every student memorizes: x = (-b ± sqrt(b^2 - 4 a c)) / (2 a). On paper it is exact. On a computer it can silently hand back a root with almost no correct digits — and the failing case is not exotic, it is whenever b^2 is much larger than 4 a c, so that sqrt(b^2 - 4 a c) is very close to |b|.
Look at the two numerators. One of them is -b + sqrt(b^2 - 4 a c) and the other is -b - sqrt(b^2 - 4 a c). When sqrt(b^2 - 4 a c) ≈ |b|, exactly one of these is a sum of two nearly-equal-magnitude numbers with opposite signs — a subtraction of near-equals — while the other is a safe sum of two same-sign numbers. The safe one is fine; the other one cancels and loses most of its digits. Crucially, it is the smaller-magnitude root that gets wrecked, which is easy to miss because the larger root still looks perfect.
Solve x^2 + 200 x + 0.000015 = 0 (a=1, b=200, c=0.000015) Exact roots: x_1 = -199.99999992500... x_2 = -7.5000000028e-8 Naive small root: x = (-b + sqrt(b^2 - 4ac)) / (2a) sqrt(40000 - 0.00006) = 199.99999985 -200 + 199.99999985 -> -0.00000015 (two large near-equals subtract!) x ~ -7.5e-8 but only ~1-2 correct digits survive Stable small root (use the good root + product of roots): x_large = (-b - sqrt(b^2 - 4ac)) / (2a) (safe: same-sign sum) x_small = c / (a * x_large) (x_1 * x_2 = c/a, no subtraction) x_small ~ -7.5000000028e-8 (full accuracy)
The one-line fix and why it works
The cure is a small algebraic rearrangement that never subtracts near-equals. First compute the root whose numerator is a safe same-sign sum — choose the sign so that you add -b and the square root with matching signs, i.e. q = -(b + sign(b)·sqrt(b^2 - 4 a c)) / 2. That gives one root reliably as x_1 = q / a. Then get the dangerous root not by the cancelling formula but from the product of roots, which for a x^2 + b x + c = 0 satisfies x_1 · x_2 = c / a. So x_2 = c / q. A division, no subtraction — the cancellation simply never happens.
- Compute the discriminant d = b^2 - 4 a c. (If d < 0 the roots are complex; if d is tiny the roots nearly coincide and that is a separate, genuine ill-conditioning of the problem, not something the rewrite can fix.)
- Form q = -(b + sign(b)·sqrt(d)) / 2. By choosing the sign to match b, you add two same-sign quantities — a safe operation with no cancellation.
- Return the two roots as x_1 = q / a and x_2 = c / q. Both are computed with only multiply and divide; neither path ever subtracts two near-equal numbers.
This is what numerical stability looks like in practice: same mathematical problem, same inputs, but a route through the arithmetic that does not amplify the unavoidable input rounding. Both the naive and the stable formula solve the identical, well-conditioned problem; only the stable one keeps the full ~16 digits, because it has a small backward error — the computed roots are the exact roots of a quadratic whose coefficients differ from yours by a hair. That is the realistic goal we aim for: not zero error, but the exact answer to a nearby problem.
The same trick everywhere else
Once you have the eye for it, you spot cancellation traps all over numerical work, and the same rewrite-the-algebra reflex defuses them. A few you will meet again on this ladder: computing 1 - cos(x) for small x (rewrite as 2 sin(x/2)^2, or use the dedicated function); the variance formula E[x^2] - E[x]^2 (use the two-pass or Welford form instead); evaluating exp(x) - 1 for small x (use expm1); and the difference quotient (f(x+h) - f(x)) / h with tiny h, where the numerator subtracts two near-equal values and the result hits a round-off floor — which is exactly why shrinking h does not keep improving a finite-difference derivative.
Two honest cautions to carry forward. First, the rewrite cures instability, never bad conditioning: if the underlying problem really is sensitive — a quadratic with a true double root, a matrix with a large condition number near 10^8 — then no algebraic cleverness recovers digits the problem itself cannot pin down; expect to lose roughly 8 of your 16. Second, every floating-point answer remains an approximation governed by machine epsilon; the goal of a stable algorithm is not exactness but to add as little as possible to the error the problem already forces on you.