Prev
Unit 1 · Computer Arithmetic and Error
Rounding Error in Arithmetic
60m100 XP
Next
🎯Computer Arithmetic and Error

Rounding Error in Arithmetic

Error of one operation, accumulation, catastrophic cancellation and compensated summation

Vancouver, January 1982
The Vancouver Stock Exchange launches a new index of its listed stocks, starting at 1000.000. It is recalculated after every trade — thousands of times a day — and each time the result is stored with three decimal places. Not rounded to three decimals. Truncated. The fourth decimal and everything after it is simply dropped. Each truncation loses less than 0.001 points; on average about 0.0005. Nobody looks at the fourth decimal of a stock index. By November 1983, twenty-two months later, the index stands at 524.811. The market has not crashed. When consultants recompute the index from the same trades with correct rounding, the answer is 1098.892. Half the value of the index was lost one ten-thousandth of a point at a time. No single operation was badly wrong.
The last lesson modelled storage: fl(x)=x(1+δ)\mathrm{fl}(x) = x(1+\delta). IEEE 754 makes the same promise for every basic operation. For ∘∈{+,−,×,÷}\circ \in \{+, -, \times, \div\} and also  \sqrt{\ }, the computed result is the exact result, correctly rounded:
fl(x∘y)=(x∘y)(1+δ),∣δ∣≤u\mathrm{fl}(x \circ y) = (x \circ y)(1 + \delta), \qquad |\delta| \le u

The standard model of floating-point arithmetic. Each operation adds at most one rounding error of relative size u — not more.

One operation is therefore never the problem. Two things are.
Accumulation. A computation with nn operations can carry nn rounding errors. If they have random signs they partly cancel, and the total typically grows like n u\sqrt{n}\,u. If they all have the same sign — truncation, or a biased algorithm — they grow like n un\,u. The Vancouver index did about a million truncations; a million times 0.0005 is about 500 points.
Cancellation. The model bounds the error of a subtraction relative to its result. When xx and yy are close, x−yx - y is tiny, and the errors xx and yy already carried — from storage, measurement or earlier operations — are now huge relative to the difference.
Derivation
How rounding errors pile up in a sum
Add three numbers in the order (a+b)+c(a + b) + c, applying the standard model to each addition, with its own rounding error δ1\delta_1 or δ2\delta_2:
s=((a+b)(1+δ1)+c)(1+δ2).s = \big((a + b)(1 + \delta_1) + c\big)(1 + \delta_2).
Multiply out and drop the products δ1δ2\delta_1\delta_2, which are of size u2u^2 and far too small to matter:
s≈(a+b+c)+(a+b)(δ1+δ2)+c δ2.s \approx (a + b + c) + (a + b)(\delta_1 + \delta_2) + c\,\delta_2.
The first numbers added carry two rounding errors, the last only one. With nn terms the pattern continues: x1x_1 and x2x_2 pass through n−1n - 1 roundings, x3x_3 through n−2n - 2, and so on, so
∣s−∑xi∣≲(n−1) u∑∣xi∣.|s - \textstyle\sum x_i| \lesssim (n - 1)\,u \sum |x_i|.
That is the worst case, when every δ\delta has the same sign — as truncation guarantees. Rounding to nearest gives δ\deltas of both signs, roughly independent of each other. Independent errors add like a random walk: their variances add, so nn errors of typical size uu give a total of typical size n u\sqrt{n}\,u. A million operations: about 106u10^6 u for a biased rule, 103u10^3 u for an unbiased one. That factor of a thousand is the whole difference between the Vancouver index and a correctly rounded one.
💡Catastrophic cancellation
Subtracting two numbers that agree in their first kk significant digits loses about kk digits of relative accuracy. The subtraction itself is exact (Sterbenz's lemma: if y/2≤x≤2yy/2 \le x \le 2y, then x−yx - y is computed without error). What it does is expose the errors already sitting in xx and yy. If each was good to 16 digits and they agree in the first 12, the difference has at most 4 good digits. The cure is never more precision first. It is algebra: rewrite the formula so the subtraction disappears.
Derivation
Why agreeing digits are lost digits
Suppose the stored values already carry small relative errors: x^=x(1+α)\hat x = x(1 + \alpha) and y^=y(1+β)\hat y = y(1 + \beta), with ∣α∣,∣β∣≤ϵ|\alpha|, |\beta| \le \epsilon. Subtract exactly:
x^−y^=(x−y)+(xα−yβ).\hat x - \hat y = (x - y) + (x\alpha - y\beta).
The second term is the error. Relative to the true difference it is
∣xα−yβ∣∣x−y∣  ≤  ϵ ∣x∣+∣y∣∣x−y∣.\frac{|x\alpha - y\beta|}{|x - y|} \;\le\; \epsilon\,\frac{|x| + |y|}{|x - y|}.
Nothing in this bound depends on how the subtraction was done; the amplification factor (∣x∣+∣y∣)/∣x−y∣(|x| + |y|)/|x - y| belongs to the numbers. If xx and yy agree in their first kk significant digits, then ∣x−y∣≈10−k∣x∣|x - y| \approx 10^{-k}|x| and the factor is about 2×10k2 \times 10^{k}: the relative error grows from ϵ\epsilon to about 10kϵ10^{k}\epsilon — exactly kk digits lost.
📐The quadratic formula, done right
For x2−108x+1=0x^2 - 10^8 x + 1 = 0 the roots are about 10810^8 and 10−810^{-8}. The textbook formula for the small root is (108−1016−4)/2(10^8 - \sqrt{10^{16} - 4})/2. In double precision 1016−4\sqrt{10^{16}-4} agrees with 10810^8 in its first 16 digits — the subtraction leaves only rounding noise, and the computed root is 7.45×10−97.45 \times 10^{-9}: 25% wrong. The stable route: compute the large root without cancellation, q=−12 (b+sign⁡(b)b2−4ac)q = -\tfrac12\,(b + \operatorname{sign}(b)\sqrt{b^2 - 4ac}), then use the product of the roots, x1x2=c/ax_1 x_2 = c/a: xlarge=q/a,xsmall=c/qx_{\text{large}} = q/a, \qquad x_{\text{small}} = c/q Now the small root is 1.0000000000000000×10−81.0000000000000000 \times 10^{-8}, correct to full precision. Same equation, same arithmetic, a different order of operations.
Derivation
Where 25% goes, and why c/qc/q is safe
The textbook route. Here a=1a = 1, b=−108b = -10^8, c=1c = 1, and b2−4ac=1081−4×10−16≈108−2×10−8\sqrt{b^2 - 4ac} = 10^8\sqrt{1 - 4 \times 10^{-16}} \approx 10^8 - 2 \times 10^{-8}. Near 10810^8 (between 2262^{26} and 2272^{27}) doubles are spaced 226×2−52=2−26≈1.49×10−82^{26} \times 2^{-52} = 2^{-26} \approx 1.49 \times 10^{-8} apart, so the square root is stored as 108−1.49×10−810^8 - 1.49 \times 10^{-8}, the nearest double. Then
xsmall=108−(108−1.49×10−8)2=7.45×10−9.x_{\text{small}} = \frac{10^8 - (10^8 - 1.49 \times 10^{-8})}{2} = 7.45 \times 10^{-9}.
The subtraction itself was exact; the answer is simply one rounding gap of the square root, halved. Every digit of the true root, 10−810^{-8}, was already lost when  \sqrt{\ } was rounded.
The stable route. If x1x_1 and x2x_2 are the roots, ax2+bx+c=a(x−x1)(x−x2)=ax2−a(x1+x2)x+ax1x2ax^2 + bx + c = a(x - x_1)(x - x_2) = a x^2 - a(x_1 + x_2)x + a x_1x_2. Comparing the constant terms: x1x2=c/ax_1 x_2 = c/a. The quadratic formula's two roots are (−b±b2−4ac)/2a(-b \pm \sqrt{b^2 - 4ac})/2a; choosing the sign that makes −b-b and the square root add (both positive here) gives the large root with no cancellation:
xlarge=−(b+sign⁡(b)b2−4ac)2a=qa.x_{\text{large}} = \frac{-\big(b + \operatorname{sign}(b)\sqrt{b^2 - 4ac}\big)}{2a} = \frac{q}{a}.
The small root then follows from the product, xsmall=ca xlarge=cqx_{\text{small}} = \dfrac{c}{a\,x_{\text{large}}} = \dfrac{c}{q} — one division, correctly rounded.
Floating-point arithmetic is considered an esoteric subject by many people.
– David Goldberg, What Every Computer Scientist Should Know About Floating-Point Arithmetic (1991)
Summation is not as simple as it looks. Adding nn numbers one at a time gives an error bound proportional to n un\,u times the size of the running sums. Once the total is large, each small new term is rounded to the total's coarse spacing and partly lost. Adding 0.1 a million times in single precision gives 100958.34 instead of 100000.
Three fixes, in increasing order of effort:
  1. Pairwise summation — add in a balanced tree. The error grows like ulog⁡2nu \log_2 n instead of u nu\,n. NumPy's np.sum does this by default.
  2. Kahan compensated summation (William Kahan, 1965) — keep a second variable that captures the low-order bits each addition drops, and feed them back into the next addition. The error becomes independent of nn to first order.
  3. Exactly rounded summation — Python's math.fsum tracks partial sums exactly and rounds once at the end.
y=xi−c,t=s+y,c=(t−s)−y,s=ty = x_i - c,\qquad t = s + y,\qquad c = (t - s) - y,\qquad s = t

One step of Kahan summation. In exact arithmetic c is always 0; in floating point (t − s) − y recovers exactly what the addition s + y rounded away.

⚠Relative and absolute error
Absolute error ∣x^−x∣|\hat x - x| has units; relative error ∣x^−x∣/∣x∣|\hat x - x|/|x| does not and counts significant digits. A position error of 1 m is excellent for a GNSS receiver on Earth's surface and useless for a docking spacecraft. Always ask which one your requirement is written in — and remember that relative error is undefined when the true value is zero.
💡Where you will use this
GNSS Navigation forms differences of nearly equal pseudoranges and carrier phases all the time; double-differencing is cancellation managed on purpose. Orbital Mechanics: the energy v2/2−μ/rv^2/2 - \mu/r is a difference of two large numbers, and its error controls orbit-determination accuracy. Structures and CFD: residual norms are long sums, where pairwise or compensated summation decides whether a solver can converge below 10−1010^{-10}.
📚Further reading
  • Nicholas Higham, Accuracy and Stability of Numerical Algorithms (2nd ed., SIAM, 2002), chapter 4, 'Summation' — every summation method in this lesson, with error bounds.
  • William Kahan, 'Further remarks on reducing truncation errors', Communications of the ACM 8(1), 40 (1965) — the one-page note that introduced compensated summation.
  • Michael Heath, Scientific Computing: An Introductory Survey (2nd ed., 2002; SIAM Classics 2018), chapter 1, 'Scientific Computing' — sources of error, cancellation and the quadratic formula.
↓But there is a problem
Rounding error comes from the machine, and we now know how to keep it small: avoid cancellation, sum carefully, prefer formulas that never subtract nearly equal numbers. But most of the errors in a numerical method are not rounding errors at all. They are put there deliberately. To compute sin⁡x\sin x the machine replaces an infinite series with a finite one. To compute a derivative it replaces a limit with a small step. In 1972 Hewlett-Packard put sine, cosine, logarithms and exponentials into a pocket calculator whose chip could do little more than add and shift. Next lesson: how — and how much error that costs.