Round-off error and catastrophic cancellation

The formula

fl(x ∘ y) = (x ∘ y)·(1 + δ), |δ| ≤ u Each single arithmetic operation (∘ is +, −, ×, ÷) lands within one unit roundoff u of the exact result. For IEEE double precision u = 2⁻⁵³ ≈ 1.1·10⁻¹⁶ (about 16 decimal digits); for single precision u = 2⁻²⁴ ≈ 6.0·10⁻⁸ (about 7 digits). One operation is essentially harmless.

relative error of (x̂ − ŷ) ≈ (x·δₓ − y·δᵧ) / (x − y) If x̂ and ŷ each arrive carrying small relative errors δₓ and δᵧ, the subtraction divides those errors by x − y. As the two values approach each other, the denominator goes to zero and the relative error of the difference goes through the roof. That is catastrophic cancellation in one line.

digits lost ≈ log₁₀( |x| / |x − y| ) A working rule of thumb: if two values agree in their first k significant digits, subtracting them throws away about k digits of precision. Agree in 8 digits under double precision and you have roughly 8 left. Agree in 6 digits under single precision and you have about one. This is the number to run in your head during a code review, the same way you'd eyeball a resistor ratio: how many leading digits do these two operands share, and how many did the format give me to start with?

Where you meet it

The DAQ post-processor that reports zero variance. The one-pass shortcut var = (Σx² − n·mean²)/(n−1) subtracts two huge, nearly identical sums. Feed it ten single-precision readings near 10,000 counts that differ in the second decimal place and it returns exactly 0.0; the true variance is about 9.2·10⁻⁴. The channel isn't dead. The formula is.
Intercept and trajectory code using the quadratic formula. When b² ≫ 4ac, whichever ± choice subtracts the square root from |b| pits two nearly equal numbers against each other. For x² − 10⁸x + 1 = 0 that's the small root (−b − √(b² − 4ac))/(2a) — here −b = +10⁸, so the minus branch is the one that cancels — and the naive formula returns 7.45·10⁻⁹ in double precision; the true value is 1.00·10⁻⁸. That's a 25 percent error from a formula everyone trusts.
Small-angle terms in pointing and alignment budgets. Compute 1 − cos(θ) for θ = 10⁻⁸ rad in double precision and you get exactly zero; the true value is 5·10⁻¹⁷. Any error budget dividing by that term just divided by zero.
Long-running fixed-point clocks. The Patriot battery at Dhahran counted time in 0.1 s ticks stored in a 24-bit register. One tenth has no exact binary representation; each tick was short by about 9.5·10⁻⁸ s. After 100 hours of continuous operation the clock had drifted 0.34 s — at Scud closing speed, a tracking-gate error of over 500 meters ^[7]^[8].

How it works

Binary floating point keeps a fixed number of significand bits, so most real numbers get rounded on the way in — including innocent-looking decimals like 0.1, which is a repeating fraction in binary. The standard model above says each operation then rounds again. Sixteen digits of headroom means individual roundings almost never matter on the bench.

Cancellation is the exception, and the mechanism is worth getting exactly right: the subtraction itself is usually exact. In IEEE arithmetic, if two floats are within a factor of two of each other, their difference is computed with zero error. Subtraction doesn't create the damage — it reveals error that was already sitting in the operands, by stripping away the matching leading digits that were correct and leaving only the trailing digits that were contaminated. Take x = 0.9876543210 and y = 0.9876542109, each trusted to 10 digits. The difference is 1.101·10⁻⁷ — three significant digits, because seven of your ten matched and vanished. Same subtraction, done in double precision, gives 1.1009999989·10⁻⁷: the machine faithfully preserved garbage.

Here is the whole disease on four lines, every number reproducible in any IEEE-754 double-precision language:

1 − cos(1e-8)                    →  0.0            (true: 5.0000·10⁻¹⁷)
2·sin²(0.5e-8)                   →  5.0000·10⁻¹⁷   (same math, rearranged)
√(10¹⁵ + 1) − √(10¹⁵)            →  1.8626·10⁻⁸    (true: 1.5811·10⁻⁸ — off 18%)
1/(√(10¹⁵ + 1) + √(10¹⁵))        →  1.5811·10⁻⁸    (conjugate form, correct)

Cancellation has a slower-burning sibling: accumulation. Summing n terms sequentially can grow error like n·u in the worst case, and in narrow formats it gets ugly fast — add 0.1 to a single-precision accumulator one million times and you get 100,958.34 instead of 100,000, high by almost one percent, because each rounded partial sum feeds the next. Sixteen digits in double precision hides this on the bench; a 24-bit format, whether float or the Patriot's fixed point, does not.

That's the split David Goldberg's survey made standard: benign cancellation subtracts exactly-known quantities and costs almost nothing; catastrophic cancellation subtracts quantities that were already approximations ^[1]^[2]. Diagnosis is therefore about the inputs, not the minus sign.

The fixes, in the order to try them:

Rearrange the algebra. Multiply by the conjugate: √(x+1) − √x = 1/(√(x+1) + √x) — the subtraction disappears. Use the stable quadratic-root pairing: compute the large root first, then x₂ = c/(a·x₁). Replace 1 − cos(θ) with 2·sin²(θ/2), which returns the correct 5·10⁻¹⁷ where the direct form returns zero.
Use the library calls built for this. log1p(x), expm1(x), and hypot(x, y) exist precisely because log(1+x), eˣ − 1, and √(x² + y²) cancel or overflow in their naive forms.
Restructure accumulations. Two-pass or Welford's method for variance; Kahan (compensated) summation when adding millions of small samples to a large running total.
Sum smart when n is large. Pairwise (tree) summation grows error like log₂(n)·u instead of n·u; sorting terms smallest-first also helps, since small values stop vanishing against a large running total.
Subtract in hardware, not software. A Wheatstone bridge exists because measuring two 5 V legs to get a 2 µV strain signal is cancellation with copper. Difference the physical quantity first, then digitize the small number.

Two boundaries on the panic. First, cancellation is a property of the formula, not the problem: the small quadratic root above is perfectly well-conditioned — a tiny change in the coefficients moves it only a little — and the stable pairing computes it to full precision. When the problem itself is ill-conditioned (a nearly singular matrix, two close eigenvalues), no rearrangement saves you, and that is a condition-number conversation, not a round-off one. Second, never test floats for equality; after any nontrivial arithmetic, compare against a tolerance scaled to the magnitudes involved, or the round-off you just read about will make if (x == y) a coin flip.

The mistakes people make: reaching for more precision instead of better algebra (double precision buys you 9 more digits than single, then the same cliff); testing numerical code only with round, friendly inputs that never trigger cancellation; and blaming the subtraction operator when the real defect is that the operands were computed sloppily upstream. And know the limit of the panic: cancellation of exact quantities is fine, and some of the best numerical tricks — compensated summation among them — deliberately use a cancelling subtraction to recover the error term. The minus sign is a tool. It just has no filter.

History

Rounding error became a subject the moment automatic computers made millions of operations cheap. In 1943, before the machines even ran, Harold Hotelling estimated that errors in Gaussian elimination could grow like 4ⁿ — a bound gloomy enough to suggest direct solvers were hopeless for large systems ^[4]^[5]. John von Neumann and Herman Goldstine answered in 1947 with a ninety-page analysis of matrix inversion in the Bulletin of the American Mathematical Society, generally counted as the founding paper of modern rounding-error analysis ^[4]^[6]. Alan Turing, who had compared notes with von Neumann at Princeton, published his own treatment in 1948 and introduced the idea of a matrix condition number along the way ^[4]^[5].

The subject's central figure was Turing's former assistant at the National Physical Laboratory, James Wilkinson, who had written floating-point subroutines for the ACE before the hardware existed. Through the 1950s and 60s Wilkinson built backward error analysis — ask not "how wrong is my answer?" but "what nearby problem did I solve exactly?" — and codified it in his 1963 book Rounding Errors in Algebraic Processes. He received the Turing Award in 1970 ^[3]^[9].

William Kahan's work drove the IEEE 754 standard of 1985, which ended the era of every vendor rounding differently and earned Kahan the 1989 Turing Award ^[10]^[11]. David Goldberg's 1991 ACM survey then gave working programmers the vocabulary — including the benign-versus-catastrophic distinction — that this entry uses ^[1]^[2]. History supplied the cautionary tale the same year: on February 25, 1991, accumulated round-off in a Patriot battery's 24-bit clock at Dhahran let a Scud through; 28 soldiers died. The GAO's report reads like a numerical analysis homework problem with casualties ^[7]^[8].

Related tools

/tools/adc-resolution — quantization is rounding; the LSB is your hardware's unit roundoff
/tools/snr-enob — effective number of bits is the measurement-chain version of "digits you actually have"
/tools/strain-gauge-bridge — microvolt signals riding on volts: cancellation handled in hardware
/tools/wheatstone-bridge — the circuit that subtracts nearly equal quantities so your code doesn't have to
/tools/vswr-return-loss — VSWR blows up as (1+|Γ|)/(1−|Γ|) when |Γ| nears 1, a live denominator-cancellation example

Sources

Written by HE in our own words from the cited sources — engineering judgment included, your stamp still required. All entries →