A real number such as \(\pi\) has infinitely many digits. Memory is finite. Something has to give. The question is what to give up so that the arithmetic is fast, predictable, and good enough for numbers as different as the mass of an electron (\(9.1\times10^{-31}\) kg) and the number of atoms in the observable universe (about \(10^{80}\)).

You already know the answer from science class. Instead of writing 0.000000000000000000000000000000911, you write \(9.11\times10^{-31}\): a handful of significant digits and an exponent that says where the decimal point goes. The decimal point "floats" to wherever it is needed; hence floating point.

A computer does the same in base 2, with a fixed budget: in the standard 64-bit "double precision" format, 53 significant binary digits (bits) and an 11-bit exponent. Fixing the number of significant digits rather than the number of decimal places has a crucial consequence: the relative accuracy is the same everywhere, but the absolute spacing between representable numbers grows with their size. Between 1 and 2 there are \(2^{52}\) evenly spaced doubles; between 1024 and 2048 there are also \(2^{52}\), spread over a range 1024 times wider.

Since 1985 essentially every processor has implemented the IEEE 754 standard. A double occupies 64 bits split into three fields.

For an ordinary ("normal") number the value is

where \(s\) is the sign bit, \(f\) is the 52-bit fraction read as an integer, and \(E\) is the 11-bit exponent field read as an integer between 1 and 2046. The leading 1 of the significand is not stored (it is always 1 for a normal number), which buys an extra bit for free: 52 stored bits give 53 bits of precision. The bias of 1023 lets the stored field be unsigned while the true exponent \(E-1023\) runs from \(-1022\) to \(+1023\).

The two leftover exponent patterns are reserved:

Let us convert 0.1 to binary by hand. To find binary digits after the point, repeatedly multiply by 2 and peel off the integer part:

So \(0.1_{10} = 0.0\overline{0011}_2 = 0.000110011001100\ldots_2\). Like \(1/3=0.333\ldots\) in decimal, it never terminates. The reason is that \(0.1 = 1/10\) and 10 has the prime factor 5, which is not a factor of the base 2. Only fractions whose denominators are powers of 2 terminate in binary.

Normalise by moving the point four places right:

In hexadecimal the whole 64-bit word is 3FB999999999999A; the final A (rather than 9) is the round-up. Because we rounded up, the stored value is slightly larger than 0.1. Python can show its exact decimal value:

Now the puzzle of the opening dissolves. The doubles nearest to 0.1, 0.2 and 0.3 are

Their exact sum is \(0.3000000000000000166533\ldots\); rounding that to the nearest double gives \(0.3000000000000000444089\ldots\), which is the double just above \(\mathrm{fl}(0.3)\). Python prints the shortest decimal string that uniquely identifies it: 0.30000000000000004.

The single most useful model of floating-point arithmetic follows from these definitions. IEEE 754 requires that the basic operations \(+,-,\times,\div\) and \(\sqrt{\ }\) be correctly rounded: the result is the exact answer rounded to the nearest double. So for any of these operations, written \(\circ\),

Read it aloud: each individual operation is wrong by a relative amount no larger than about \(10^{-16}\). All the drama in this chapter comes not from a single operation but from how such tiny relative errors are amplified by what we do next.

The diagram shows the spacing between doubles as a function of their size. It is a staircase (constant within each power-of-two band, then doubling) that looks like a straight line on log-log axes. At \(2^{53}=9\,007\,199\,254\,740\,992\) the gap reaches 1, so \(2^{53}+1\) cannot be stored: float(2**53 + 1) returns 9007199254740992.0. This is why JavaScript, whose only number type is the double, defines Number.MAX_SAFE_INTEGER as \(2^{53}-1\).

Solve \(x^2 + 10^8x + 1 = 0\). The quadratic formula from Chapter 8 gives Try it: Cancellation and the stable rewrite.

The two roots are about \(-10^8\) and about \(-10^{-8}\) (their product is \(c/a=1\)). In double precision:

The large root is fine. The small root should be \(-1.0000000000000000\ldots\times10^{-8}\); we got \(-7.45\times10^{-9}\), wrong in the first significant digit. A 25% error from a calculation in which every single operation was accurate to 16 digits. What happened?

Here \(b=10^8\) and \(\sqrt{b^2-4}\) are equal in their first 16 digits:

The computed square root carries a relative error up to \(u\), i.e. an absolute error up to about \(10^8\times10^{-16}=10^{-8}\). Subtraction \(-b+d\) is itself performed exactly (the two numbers are so close that the difference is representable), but it removes the 16 digits they share and leaves only what was in the "noise" region. The absolute error of about \(10^{-8}\) is now attached to a result of size \(2\times10^{-8}\): relative error of order 1. This is catastrophic cancellation.

The key insight: subtraction does not create error; it exposes error that was already present in its inputs. If the inputs were exact, the difference would be exact.

The cure is algebra. Rationalise the numerator (multiply top and bottom by the conjugate):

Now the denominator adds two numbers of the same sign: no cancellation. With our numbers, \(2c/(-b-d) = 2/(-2\times10^8) = -1.0\times10^{-8}\) to full precision. Equivalently, use the product of roots: \(x_{\text{small}} = c/(a\,x_{\text{large}})\).

Some standard rewrites:

For \(x=10^{-8}\), 1 - math.cos(x) returns exactly 0.0 (because \(\cos(10^{-8}) = 1 - 5\times10^{-17}\) rounds to 1), while 2*math.sin(x/2)**2 returns 5.0000000000000005e-17, correct to 16 digits.

The other two floating-point hazards are about range, not precision.

Overflow. Computing \(\sqrt{x^2+y^2}\) for \(x=10^{200}\) overflows at \(x^2\) even though the answer, about \(10^{200}\), is perfectly representable. math.hypot rescales internally.

Underflow. Multiply the probabilities of 1000 independent events, each 0.01, and the true answer \(10^{-2000}\) underflows to 0. Every probabilistic model that multiplies many likelihoods therefore works with logarithms: \(\log\prod p_i = \sum\log p_i\) turns a product that underflows into a sum that is perfectly well behaved (Chapter 12).

When a computed answer is bad, whose fault is it? Sometimes the problem is so sensitive that no algorithm could do better; sometimes a perfectly reasonable problem was solved with a bad algorithm. Numerical analysts keep these two ideas strictly separate, and you should too.

Consider the system Try it: Nudge the right-hand side; Feel κ ≈ 4 × 10⁴.

Subtracting the first equation from the second gives \(0.0001\,y=0.0001\), so \(y=1\) and \(x=1\). Now change the right-hand side of the second equation by one part in 20,000, from \(2.0001\) to \(2.0002\). Then \(0.0001\,y = 0.0002\), so \(y=2\) and \(x=0\). A change in the fifth significant digit of the data moved the answer from \((1,1)\) to \((0,2)\).

Geometrically, the two equations are lines that are almost parallel (Chapter 6). Nudge either one and their intersection slides a long way along them. No algorithm, however careful, can fix that: if the data carry an uncertainty of \(10^{-4}\), the answer is uncertain in its first digit.

Rule of thumb: you can lose about \(\log_{10}\kappa\) significant digits. With 16 digits available and \(\kappa=4\times10^4\), expect about 11 or 12 good digits from any decent solver, if the data themselves are exact.

The function \(f(x)=x-1\) near \(x=1\) is a one-variable example: \(\kappa = |x\cdot1/(x-1)|\), which blows up as \(x\to1\). That is catastrophic cancellation seen through the lens of conditioning: subtracting nearly equal numbers is an ill-conditioned operation.

An algorithm is stable if it does not add much error beyond what the conditioning of the problem forces. The gold standard is backward stability: the computed answer is the exact answer to a slightly perturbed problem (perturbed by a few multiples of \(u\)). Then

This one line is the most important idea in numerical analysis. It splits every error into "the problem’s sensitivity" times "the algorithm’s sloppiness".

Define \(I_n = \int_0^1 x^n e^{x-1}\,dx\) for \(n=0,1,2,\ldots\). Integration by parts (Chapter 19) with \(u=x^n\), \(dv=e^{x-1}dx\) gives a recurrence: Try it: Forwards versus backwards recurrence.

Every \(I_n\) is positive (the integrand is positive) and they decrease (\(x^n\) shrinks on \(\lbrack 0,1\rbrack\)). Run the recurrence forwards in Python:

A negative value at \(n=18\), then nonsense. Why? \(I_0\) is stored with error \(\delta\approx10^{-17}\). Each step multiplies the inherited error by \(-n\), so after 20 steps it has been multiplied by \(20!\approx2.4\times10^{18}\). The algorithm is wildly unstable, even though the problem (compute \(I_{20}\)) is perfectly well-conditioned.

Now run the same recurrence backwards: \(I_{n-1} = (1-I_n)/n\). Each step divides the error by \(n\). Start at \(n=40\) with the crude guess \(I_{40}=0\) (true value about 0.024, so the starting error is about 0.024):

The starting error has been divided by \(40\cdot39\cdots21\approx 3.4\times10^{29}\); the answer is correct to all 16 digits. Same mathematics, same recurrence, opposite direction, opposite fate.

Solve \(f(x)=0\) when algebra gives up: \(\cos x = x\), \(x^5 - x - 1 = 0\), or "what interest rate makes this loan’s payments equal its present value?" You met Newton’s method in Chapter 18. Here we put it alongside its rivals and compare them as algorithms: what they need, how fast they converge, and how they fail.

To measure errors exactly we use a test problem whose answer we know: \(f(x)=x^2-2\), root \(r=\sqrt2 = 1.41421356237309\ldots\).

Intuition. If a continuous function is negative at \(a\) and positive at \(b\), it must cross zero in between (the Intermediate Value Theorem, Chapter 16). Check the midpoint; keep whichever half still has a sign change. It is binary search on the real line. Try it: Race all four methods.

By hand on \(\lbrack 1,2\rbrack\) (\(f(1)=-1<0\), \(f(2)=2>0\)):

Each step halves the bracket, so after \(n\) steps the error is at most \((b-a)/2^{n}\). Reaching \(10^{-10}\) from a bracket of width 1 needs \(n\ge\log_2(10^{10})\approx33.2\), i.e. 34 steps. Linear convergence with ratio \(\tfrac12\): about one decimal digit every 3.3 steps.

Newton replaces the curve by its tangent line at \(x_n\) and jumps to where the tangent crosses zero:

For \(f(x)=x^2-2\) this is \(x_{n+1}=x_n-\dfrac{x_n^2-2}{2x_n} = \dfrac{x_n}{2}+\dfrac{1}{x_n}\). From \(x_0=1\):

The errors are \(8.6\times10^{-2}\), \(2.5\times10^{-3}\), \(2.1\times10^{-6}\), \(1.6\times10^{-12}\): the exponent roughly doubles each time. That is quadratic convergence. Step 5 hits the correctly rounded value of \(\sqrt2\). (Step 6 can bounce to the neighbouring double and back; at the level of the last bit, rounding decides, not Newton. Stopping tests must allow for that.)

Why quadratic? Expand \(f\) about the root with Taylor’s theorem: the tangent line captures the first-order term exactly, so the leftover error comes from the second-order term, which is proportional to \(e_n^2\). Precisely, \(e_{n+1}\approx \left|\dfrac{f''(r)}{2f'(r)}\right| e_n^2\).

Often \(f'\) is unavailable (\(f\) is a simulation, or a black box). Replace the tangent by the secant through the last two points, i.e. approximate \(f'(x_n)\) by \(\dfrac{f(x_n)-f(x_{n-1})}{x_n-x_{n-1}}\): Try it: Step the secant method.

From \(x_0=1, x_1=2\) (so \(f(x_0)=-1, f(x_1)=2\)):

Errors: \(8.1\times10^{-2}, 1.4\times10^{-2}, 4.2\times10^{-4}, 2.1\times10^{-6}, 3.2\times10^{-10}\). Each exponent is about 1.6 times the previous one: the order is the golden ratio \(p=\frac{1+\sqrt5}{2}\approx1.618\). Slower per step than Newton, but each step costs one new function evaluation instead of a function plus a derivative, so per unit of work the secant method is often faster.

Rewrite \(f(x)=0\) as \(x=g(x)\) and iterate \(x_{n+1}=g(x_n)\). A solution of \(x=g(x)\) is a fixed point of \(g\). Newton’s method is itself fixed-point iteration with \(g(x)=x-f(x)/f'(x)\). Try it: Watch x ← 2/x oscillate; Fixed point with cobweb.

Why it works. Near the fixed point \(r\), \(g(x_n)-g(r)\approx g'(r)(x_n-r)\), so \(e_{n+1}\approx|g'(r)|\,e_n\). The iteration converges (for starting points close enough) when \(|g'(r)|<1\), linearly with ratio \(|g'(r)|\), and the smaller \(|g'(r)|\), the faster.

Two rearrangements of \(x^2=2\):

Fixed-point thinking reappears everywhere: Jacobi iteration below, PageRank (Chapter 28), value iteration in reinforcement learning, and the convergence proofs of many optimisers.

Decision rule. If you can bracket the root and you only need it once, bisection’s certainty is worth its slowness. If you need the root to full precision, or millions of times, use a fast method, ideally safeguarded: production libraries (for example scipy.optimize.brentq) use Brent’s method, which tries secant-like and inverse-quadratic steps but falls back to bisection whenever a fast step would leave the bracket. You get guaranteed convergence and superlinear speed. Use Newton when derivatives are cheap and a good starting point is known (square roots in hardware, inside optimisers, implicit ODE solvers).

Here are all four, ready for experiments:

Bisection and fixed-point iteration stop once their steps fall below \(10^{-12}\), so they agree with \(\sqrt2\) to about 12 digits, exactly as their tolerances promise. Secant reaches the correctly rounded value; Newton finishes one ulp below it (the last-bit bounce mentioned earlier). Each method delivers what its stopping test asked for, and no more.

You have a function as code, say a physics simulation or a pricing model, and need its derivative. The definition (Chapter 17) suggests taking a small \(h\):

Mathematics says "the smaller \(h\) the better". The computer disagrees.

Truncation error comes from using a finite \(h\) instead of a limit. Taylor’s theorem (the polynomial approximation of Chapter 18) says

so

Truncation error is proportional to \(h\): it shrinks as \(h\) shrinks.

Rounding error comes from the subtraction \(f(x+h)-f(x)\): two nearly equal numbers, each carrying a rounding error of about \(u|f(x)|\). That is catastrophic cancellation again. The absolute error of the numerator is about \(u|f(x)|\), and we then divide by \(h\), giving an error of about \(u|f(x)|/h\), which grows as \(h\) shrinks.

Total error is approximately

Minimise with calculus: \(E'(h)=\tfrac{M}{2}-\tfrac{u|f|}{h^2}=0\) gives \(h^\ast=\sqrt{2u|f|/M}\), which is of order \(\sqrt{u}\approx10^{-8}\) when \(f\) and \(f''\) are of order 1. At that \(h\) both terms are about \(10^{-8}\): the forward difference can only ever deliver about half the available digits.

The central difference is smarter. Subtract the Taylor expansions of \(f(x+h)\) and \(f(x-h)\); the \(f''\) terms cancel:

Truncation error is now \(O(h^2)\). Balancing \(h^2\) against \(u/h\) gives \(h^\ast\sim u^{1/3}\approx 6\times10^{-6}\) and a best error around \(u^{2/3}\approx10^{-11}\): about two-thirds of the digits.

Take \(f(x)=\sin x\) at \(x=1\), where the true derivative is \(\cos1 = 0.5403023058681398\). Try it: Plot the V-shaped error.

Read the forward column top to bottom: the error falls by a factor of 10 each time \(h\) does (order 1), bottoms out near \(h=10^{-8}\), then rises again. The central column falls by a factor of 100 per decade of \(h\) (order 2), bottoms out near \(h=10^{-6}\), then rises. Plot the error against \(h\) on log–log axes and you get the characteristic V.

Many integrals have no elementary antiderivative: \(\int e^{-x^2}dx\) (the normal distribution, Chapter 30), \(\int\sqrt{1+\cos^2x}\,dx\) (arc length of a sine wave), and every integral of measured data. The Riemann sums of Chapter 19 converge, but slowly. Can we do better with the same number of function evaluations?

As a test integral with a known answer we use

Split \(\lbrack a,b\rbrack\) into \(n\) strips of width \(h=(b-a)/n\), with nodes \(x_i=a+ih\).

Trapezoid rule. Replace the curve on each strip by the straight line joining its endpoints. Each strip is a trapezoid with area \(\tfrac h2(f(x_{i})+f(x_{i+1}))\). Adding them, every interior node is counted twice:

The error is \(-\dfrac{(b-a)h^2}{12}f''(\xi)\) for some \(\xi\): order 2. (The sign tells you something too: for a convex function, \(f''>0\), the chords lie above the curve and the trapezoid rule overestimates.)

Simpson’s rule. Replace the curve on each pair of strips by the parabola through three points. Integrating that parabola (an exercise in Chapter 19-style algebra) gives weights \(\tfrac h3(1,4,1)\) on each pair, and adding pairs (\(n\) even):

The error is \(-\dfrac{(b-a)h^4}{180}f^{(4)}(\xi)\): order 4. Surprisingly, it is exact for cubics, not just quadratics, because the cubic error term is odd about each pair’s midpoint and cancels.

With \(n=2\), \(h=0.5\), nodes \(0,0.5,1\):

Trapezoid:

Error \(1.7539310925-1.7182818285=0.0356\), an overestimate, as predicted for convex \(e^x\).

Simpson:

Error \(5.8\times10^{-4}\), sixty times smaller with the same three function values.

Doubling \(n\) repeatedly confirms the orders:

Trapezoid errors drop by a factor of 4 (\(=2^2\)) per doubling; Simpson errors by 16 (\(=2^4\)).

A completely different idea from probability (Chapter 29): the average value of \(f\) on \(\lbrack a,b\rbrack\) is \(\frac{1}{b-a}\int_a^b f\), so pick \(n\) uniformly random points and average:

By the central limit theorem (Chapter 30) the error is random with standard deviation \(\sigma/\sqrt n\), where \(\sigma\) is the standard deviation of \(f(U)\). To gain one decimal digit you need 100 times as many points. Dreadful in one dimension, as the plot shows. But the \(1/\sqrt n\) rate does not depend on the dimension, whereas a grid-based rule in \(d\) dimensions with \(m\) points per axis needs \(m^d\) evaluations. For a 10-dimensional integral (pricing a basket option, rendering a pixel with light bouncing several times), Monte Carlo wins easily.

Systems \(A\mathbf x=\mathbf b\) with thousands or millions of unknowns arise in structural engineering, circuit simulation, least-squares fitting (Chapter 33), and the discretisation of every partial differential equation. You know Gaussian elimination from Chapter 6 and Chapter 27. What changes when it runs in floating point?

Solve Try it: Eliminate with a 10⁻²⁰ pivot.

The exact solution is extremely close to \(x=1, y=1\). Naive elimination uses the tiny \(10^{-20}\) as the pivot: multiply row 1 by \(m=1/10^{-20}=10^{20}\) and subtract from row 2:

\(x=0\) is completely wrong. The huge multiplier swamped row 2’s information (the "1" and "2" were rounded away), and the back substitution then divided a cancellation error by a tiny number.

Partial pivoting swaps the rows first so that the pivot is the largest entry (in absolute value) in its column:

Now the multiplier is \(10^{-20}\), at most 1 in size, and nothing gets swamped. With partial pivoting every multiplier satisfies \(|m|\le1\), which keeps the growth of rounding errors under control; in practice Gaussian elimination with partial pivoting is backward stable.

Solve

Column 1. The entries are \(2,-3,-2\); the largest in absolute value is \(-3\) (row 2). Swap rows 1 and 2:

Multipliers: \(m_{21}=2/(-3)=-\tfrac23\), \(m_{31}=-2/(-3)=\tfrac23\).

Column 2. Among rows 2 and 3 the entries are \(\tfrac13\) and \(\tfrac53\); pivot on \(\tfrac53\), so swap rows 2 and 3. Multiplier \(m_{32}=\tfrac13\big/\tfrac53=\tfrac15\):

The triangular system is

Back substitution.

Check in the original first equation: \(2(2)+3-(-1)=8\). Correct.

In practice call numpy.linalg.solve, which calls the LAPACK library’s pivoted LU factorisation. Never compute \(A^{-1}\) and multiply by it: that costs about three times as much and is less accurate.

For a sparse system with a million unknowns, each equation involving only a handful of neighbours (a grid in a heat-flow simulation, links in a web graph), elimination is hopeless: it fills in the zeros and costs \(\sim n^3\). Iterative methods instead start from a guess and improve it, touching only the nonzero entries each sweep.

Jacobi’s idea: solve equation \(i\) for unknown \(i\), using the previous values of all the others. For

Start from \((x_0,y_0)=(0,0)\):

The iterates spiral in on the exact solution \((2,1)\); by \(k=8\) they read \((1.99990, 0.99995)\). Each sweep shrinks the error by a roughly constant factor (here about \(1/\sqrt{12}\approx0.29\)): linear convergence, exactly the fixed-point iteration of the root-finding section in vector form.

When does it converge? A sufficient condition is that \(A\) is strictly diagonally dominant: in each row, the diagonal entry is larger in absolute value than the sum of the absolute values of the others (\(4>1\), \(3>1\) above). Swap the equations to \(x+3y=5,\ 4x+y=9\) and Jacobi diverges.

Gauss–Seidel uses each new value as soon as it is computed and typically converges about twice as fast; conjugate gradient (for symmetric positive-definite \(A\)) and GMRES are the modern workhorses.

Mechanics students learn Newton’s second law: force equals mass times acceleration. For a pendulum of length \(L\) swinging through angle \(\theta(t)\) this gives

For small swings, \(\sin\theta\approx\theta\) and the solution is a cosine (Chapter 13). For large swings there is no elementary formula. Planetary orbits, epidemics, circuits and chemical reactions are the same story: an equation for the rate of change of a state, and no formula for the state itself. We step forward in time instead.

Intuition. \(F(t,y)\) tells you the slope of the solution right now. Walk in that direction for a time \(h\), then look again:

It is the tangent-line approximation of Chapter 17 applied over and over.

By hand. Take \(y'=y,\ y(0)=1\), whose exact solution is \(e^t\), and step to \(t=1\) with \(h=0.25\): Try it: Try Euler with h = 0.25.

The true value is \(e = 2.71828\); Euler is low by 0.277 because the slope keeps increasing during each step and Euler always uses the stale slope from the start. (Notice that \(y_4=1.25^4\): Euler computes \((1+h)^{1/h}\), the compound-interest approximation of \(e\) from Chapter 12.) With \(h=0.1, 0.01, 0.001\) the errors are \(0.125, 0.0135, 0.00136\): the error is proportional to \(h\). Euler is a first-order method; for six correct digits you would need about a million steps.

Intuition. Instead of trusting the slope at the start of the step, probe it at the start, twice at the midpoint, and at the end, and take a weighted average, in exactly the proportions of Simpson’s rule (\(1:2:2:1\) is \(1:4:1\) with the middle weight split):

By hand, one step of \(y'=y\), \(y_0=1\), \(h=0.25\):

The exact value is \(e^{0.25}=1.2840254\). One RK4 step is already correct to 5 digits, where Euler’s first step (1.25) was correct to 1. Four steps reach \(t=1\) with error \(7.2\times10^{-5}\), compared to Euler’s 0.277. RK4 is fourth order: halve \(h\) and the error falls by 16. It costs four evaluations of \(F\) per step, but that is a bargain.

A frictionless pendulum conserves energy: per unit mass and with \(g/L=1\), Try it: Run the energy-drift experiment.

(kinetic plus potential) is constant forever. Energy therefore makes a superb diagnostic: any drift in \(E\) is pure numerical error. Release the pendulum from rest at \(\theta=1\) radian and simulate 60 time units with \(h=0.1\), using Euler, RK4, and a third contender, the semi-implicit (symplectic) Euler method, which differs from Euler only in updating the velocity first and then using the new velocity to update the position:

Try it: Run this in your browser.

Three very different personalities:

The lesson: order of accuracy is not the only virtue. For long simulations of conservative systems (planetary orbits over millions of years, molecular dynamics, game physics), preserving the right structure beats raw accuracy per step.

Add 0.1 to itself ten million times. The answer should be 1,000,000. Try it: Add 0.1 ten million times.

Not only is 0.1 inexact, but once the running sum is large each addition rounds the tiny increment to the ulp of the sum, and the rounding errors accumulate. For naive summation of \(n\) terms the error bound grows like \(n\,u\sum|x_i|\): with \(n=10^7\) that is a loss of 7 of the 16 digits in the worst case.

Imagine a cash register that can only display four significant digits. The total is 1000; a customer adds a 0.4 item; the register shows 1000 again and the 0.4 vanishes. Ten such items vanish entirely: the register says 1000 when the truth is 1004. A sensible cashier would keep the lost amounts in a jar and add them back in when they became big enough to register. That jar is Kahan’s compensation variable. Try it: Run the cash register.

Here is the cash register in Python, using the decimal module set to 4 significant digits so we can watch every step:

Trace the first few Kahan steps after the sum reaches 1000:

\(c\) always holds the negative of the rounding error just committed (in \(t-s\), the part of \(y\) that actually made it into the sum, minus \(y\)). Algebraically \((t-s)-y=0\); in floating point it is exactly the lost low-order part. This is one of the rare places where you want the subtraction of nearly equal numbers, because here they are exact.

You have values of a function only at a few points: a table of tax brackets, a few calibration measurements, keyframes of an animation. You need values in between. Interpolation finds a function that passes exactly through the given points and uses it to estimate the rest.

Through any two points there is exactly one line; through three points (with distinct \(x\)) exactly one parabola; through \(n+1\) points exactly one polynomial of degree at most \(n\) (Chapter 7). Lagrange’s idea writes it down directly. For nodes \(x_0,\ldots,x_n\), build basis polynomials that are 1 at their own node and 0 at all the others:

The \(\prod\) symbol is the product analogue of \(\sum\): multiply the factors for every \(j\) from 0 to \(n\) except \(j=i\). At \(x=x_i\) every factor equals 1; at any other node \(x_k\) the factor with \(j=k\) is zero. Then

passes through every data point, because at \(x_k\) only the \(k\)-th term survives and it equals \(y_k\).

Worked example. Find the polynomial through \((0,1),(1,3),(2,7)\) and estimate \(p(1.5)\).

Check: \(p(0)=1\), \(p(1)=3\), \(p(2)=7\). And \(p(1.5)=2.25+1.5+1=4.75\). (Directly from the basis values: \(L_0(1.5)=-0.125\), \(L_1(1.5)=0.75\), \(L_2(1.5)=0.375\), so \(p(1.5)=-0.125+2.25+2.625=4.75\).)

More points should mean a better fit. Carl Runge showed in 1901 that this is false. Interpolate \(f(x)=1/(1+25x^2)\) on \(\lbrack -1,1\rbrack\) at equally spaced points: Try it: See Runge’s wiggle; Same degree, Chebyshev nodes.

With equally spaced nodes the interpolant fits beautifully in the middle and oscillates wildly near the ends, and the oscillation grows with the degree. The cure is either to cluster the nodes near the ends, using Chebyshev nodes \(x_k=\cos\!\big(\tfrac{(2k+1)\pi}{2(n+1)}\big)\) (the projections onto the \(x\)-axis of equally spaced points on a semicircle), or to abandon a single high-degree polynomial altogether.

The order matters. There is no point computing to 15 digits a model that is right to 2, or polishing rounding error when the step size is the dominant error.