Scientific Computing · Floating Point
Catastrophic Cancellation and the Cases Where Subtraction Is Exact
Why subtracting two nearly equal numbers amplifies existing error rather than creating it, when the difference is exact by Sterbenz, and how formulas are rewritten to avoid the loss.
The textbook warning is that subtracting nearly equal numbers destroys precision. The phrasing is misleading in a way that matters for diagnosis. Floating-point subtraction is one of the most accurate operations available: when the exact difference of two representable values is itself representable, the hardware returns it with no error at all. What cancellation destroys is not precision in the subtraction but the credibility of the operands, and the arithmetic is merely the messenger.
An operation that is exact whenever it can be
Correct rounding has a consequence that is easy to miss. If the exact result of an operation is representable in the destination format, the correctly rounded result is that value — there is nothing to round. For addition and subtraction of two floating-point numbers, this means the computed result is exact whenever the exact result fits.
That is a strong statement, and it makes subtraction unusual among the four operations. A product of two p-bit significands generally needs 2p bits, so it is almost always rounded. A quotient is almost never exact. A sum or difference of two numbers of similar magnitude has a short significand and usually fits. In binary32, the difference of two values in [1, 2) is a multiple of 2^-23 and less than 1, so it needs at most 23 bits and is always exactly representable.
The general condition is Sterbenz’s lemma: if y/2 ≤ x ≤ 2y for positive x and y, then x − y is exactly representable. The proof is a counting argument. Both operands are multiples of 2^-t scaled by their common binade, so the difference is a multiple of the smaller spacing; and the difference is smaller than y, so its significand needs no more bits than y already has. The lemma extends to the subnormal range, where the spacing is uniform.
Two consequences follow. First, exact cancellation is not an error event: if x and y are exact inputs within a factor of two, x − y is exact, and a program that computes 1.5 - 1.25 and gets 0.25 has lost nothing. Second, when cancellation does hurt, the damage was already present in the operands, and the subtraction has only rescaled it.
Error that arrives with the operands
Suppose x̂ and ŷ are computed values with relative error bounded by u, so x̂ = x(1 + δ₁) and ŷ = y(1 + δ₂) with |δᵢ| ≤ u. Their difference is
$$ \hat{x} - \hat{y} = (x - y) + (x\delta_1 - y\delta_2) $$
The absolute error is bounded by u(|x| + |y|), which does not depend on how close the operands are. The relative error of the difference is therefore
$$ \frac{u(|x| + |y|)}{|x - y|} $$
and that quantity has no bound as x → y. This ratio is the condition number of subtraction, and it is the only number needed to predict the damage.
Make it concrete. Take two binary32 values near 1, each carrying a relative error of 2^-24. Their absolute error is about 6e-8 each, so up to 1.2e-7 combined. If the true difference is 2^-23, about 1.19e-7, the computed difference has a relative error near 1 — every digit is wrong. If the true difference is 0.5, the relative error is about 2.4e-7, which is fine. Same operands, same subtraction, completely different reliability, and the only difference is the magnitude of the result relative to the operands.
A useful diagnostic follows: after a subtraction, compare the exponent of the result with the exponent of the larger operand. A drop of many binades means the result’s relative accuracy is that much worse than the operands’. If the operands were themselves the output of a longer computation, that loss compounds.
Reformulation instead of higher precision
The standard remedy is algebraic, not arithmetic. If a formula subtracts two quantities that are close by construction, there is usually an identity that computes the same value without the subtraction, and it is usually cheaper than widening the precision.
The classic cases are all in the same family, where a function is evaluated near a point where two terms nearly cancel:
1 − cos xfor smallxloses accuracy quadratically. The absolute error incos xis aboutu, while the difference is aboutx²/2, so the relative error of the result is roughly2u/x²— atx = 2^-20in binary64 that is a relative error near2^-12from a single subtraction. The identity2 sin²(x/2)computes the same quantity with no cancellation. The C library exposesexpm1andlog1pfor the two most common instances,e^x − 1andlog(1 + x), precisely because the naive forms fail for small arguments.√(x+1) − √xfor largexhas a relative error that grows in proportion tox, since the difference is about1/(2√x)while the absolute error in each square root is aboutu√x. Multiplying by the conjugate gives1/(√(x+1) + √x), in which every term is of the same magnitude.- The quadratic formula loses one root when
b² ≫ 4ac. Withbpositive,−b + √(b² − 4ac)subtracts two nearly equal numbers, while the other root,−b − √(b² − 4ac), is computed accurately.
The last case is worth writing out because the fix is not an approximation but an exact rearrangement. The product of the two roots is c/a, so once one root is known accurately the other follows by division.
#include <math.h>
/* Roots of a*x^2 + b*x + c. The root that would cancel is obtained
from the product of the roots instead of from the difference. */
static void quadratic(double a, double b, double c, double *r1, double *r2)
{
double disc = b * b - 4.0 * a * c;
if (disc < 0.0) { *r1 = *r2 = NAN; return; }
double s = sqrt(disc);
double q = -0.5 * (b + copysign(s, b)); /* same sign as b: terms add */
*r1 = q / a;
*r2 = (q != 0.0) ? c / q : 0.0;
}
copysign selects the sign that makes the two terms add rather than cancel. When b is zero the choice is arbitrary and both roots are symmetric; when c is zero, q is zero and the second root is genuinely zero.
Cancellation that is not written as subtraction
Most damaging cancellation is hidden inside a formula that contains no visible minus sign between close quantities. Three patterns account for a large share of it.
The first is a sum of terms with mixed signs, which is the general case of which subtraction is the two-term instance. A dot product Σ xᵢyᵢ has condition number Σ|xᵢyᵢ| / |Σ xᵢyᵢ|; when the terms cancel heavily, the condition number is large and the accumulated rounding error is amplified by exactly that factor. The same expression governs the accuracy of any reduction with signed contributions.
The second is the variance computed from raw moments. The identity Var = E[x²] − (E[x])² subtracts two large, nearly equal quantities whenever the standard deviation is small compared with the mean, and the condition number grows as the square of that ratio. Shifting the data by a constant before squaring removes most of the problem, and the two-pass algorithm — compute the mean, then sum the squared deviations — removes the rest at the cost of a second pass over the data. Welford’s incremental update avoids the second pass while remaining stable, which is why it appears in streaming statistics libraries.
The third is evaluation of a polynomial near a root, or of any function near a zero crossing. Horner’s rule evaluates p(x) accurately in the sense of backward error, but the value near a root is a small difference of large terms, so the relative accuracy of the root location is limited by the conditioning of the polynomial, not by the evaluation. A root that is ill-conditioned cannot be recovered by evaluating more carefully; the polynomial’s coefficients have to be reparameterised.
What the diagnosis implies
Cancellation is a property of the problem, not of the arithmetic. The condition number of a subtraction is determined by the data, and no implementation strategy removes it: evaluating the same expression in binary128 postpones the failure but does not change the ratio (|x|+|y|)/|x−y|. Wider precision helps only when the error in the operands comes from rounding rather than from the data, which is a different situation — the operands are then accurate to 2^-113 instead of 2^-53, and the amplified error is correspondingly smaller.
Distinguishing the two situations is the practical skill. If the operands are exact inputs, subtraction is exact and the result is trustworthy. If the operands are computed, the error in the result is the error in the operands amplified by the condition number, and the correct response is either a reformulation or a wider format, chosen by asking where the operands’ error came from. That question is the subject of Forward Error, Backward Error, and Condition Numbers, and the summation case, where cancellation accumulates over many terms, is handled by the techniques in Compensated Summation and How Summation Order Bounds Error.
References
- Accuracy and Stability of Numerical Algorithms, 2nd edition — Nicholas J. Higham; SIAM
- Floating Point Computation — Pat H. Sterbenz; Prentice-Hall, 1974
- IEEE 754-2019: Standard for Floating-Point Arithmetic — IEEE