Where floating point bites

Nine small experiments, each a tab, that show where computer arithmetic parts company with real arithmetic — and the standard cure for each. Everything is computed live in your browser, in IEEE double precision unless a tab says otherwise; float32 is emulated with Math.fround and float16 with a half-precision rounding helper, rounding after every operation.

The tabs

Derivatives — slide log₁₀ h from 0 to −20 and watch the error of the forward difference, the central difference and the complex step for sin x, eˣ, x³ or √x. Solid curves are computed in the chosen precision; dotted curves are the same formulas in exact arithmetic, so you can see where rounding takes over. The bits of x and x + h are shown with the differing bits highlighted, until x + h == x.
Quadratic — x² + bx + 1 (or any a, b, c) solved by the textbook formula and by the stable formula q = −(b + sign(b)√D)/2, roots q/a and c/q, with each digit green if correct and red if not, in 16, 7 (float32) or 4 significant digits. The plot shows the small root’s relative error as |b| grows from 1 to 10¹⁰.
(1 + 1/n)ⁿ — approaching e, drifting away, then collapsing to exactly 1 once 1 + 1/n rounds to 1; with float32 and the log1p fix for comparison.
Sums — Kahan’s cash register on a 4-digit machine, traced step by step; 0.1 added up to ten million times (naive, Kahan, and the exact decimal answer); and 1,000 floats summed in different orders, each total measured in ulps from the correctly rounded sum.
Horner — expanded (x − 1)ⁿ evaluated near x = 1 by the naive loop, a running power and Horner’s rule, against the factored form, with multiplication counts.
Recurrence — Iₙ = 1 − n·Iₙ₋₁ run forwards (it goes negative at n = 18) and backwards from a crude guess (correct to every digit).
Conditioning — the nearly parallel lines x + y = 2, x + (1 + ε)y = 2 + ε: nudge the right side and watch the intersection slide, with the condition number κ₂ and the measured amplification.
Pivoting — the 2×2 system with a 10⁻²⁰ pivot, eliminated with and without partial pivoting.
Runge — interpolate 1/(1 + 25x²) through 2 to 25 equally spaced or Chebyshev nodes, with one polynomial or a cubic spline; drag a data point to see how far the disturbance spreads.

Limits

float16 and float32 are emulated: each elementary result (sin, exp, sqrt, +, ×, …) is computed in double and then rounded, which is what a correctly rounded float16/32 library would give. Hardware float16 maths libraries may differ in the last bit.
In the rounding-dominated part of a V curve the exact noise depends on the last bit of the platform’s sin and exp, so a few values there can differ from the book’s Python output (for example the forward difference at h = 10⁻⁸); the shape and the best h do not.
The quadratic tab’s true roots are computed exactly from the stored coefficients with 90-digit integer arithmetic. The 4-digit mode rounds every operation to 4 significant digits (round half up), like a by-hand calculator.

Related: Root finding, step by step (bisection, Newton, secant, fixed point), Graphing calculator, Linear systems, visually (row operations and nearly parallel lines), Matrix calculator.

Companion to Mathematics, Reconnected: Numerical methods (floating point, cancellation, conditioning, the V-shaped error, pivoting, Kahan summation, interpolation), with earlier glimpses in Numbers, Fractions, Systems of equations, Polynomials, Quadratics, Limits and Derivatives.

Linking to an example

Every setting can be written by hand in the address, for example /tools/numerical-error/?b=1e8&digits=16#quad. The part after # picks the tab: #diff, #quad, #limit, #sum, #horner, #recur, #cond, #pivot or #runge. Settings a link leaves out take their defaults. Numbers may be written as 1e8, -1e8, 0.001.

Parameter Meaning Example

Parameter	Meaning	Example
`fn`	derivatives: `sin`, `exp`, `cube` (x³) or `sqrt`	`sin`
`x`	derivatives: the point (default 1, or 0.01 for √x)	`1`
`prec`	derivatives: `f16`, `f32` or `f64`	`f32`
`h`	derivatives: log₁₀ of the step h, −20 to 0	`-8`
`a`, `b`, `c`	quadratic: coefficients of ax² + bx + c	`1`, `-1e8`, `1`
`digits`	quadratic: significant digits, `16`, `7` or `4`	`4`
`n`	(1 + 1/n)ⁿ: the value of n	`1e15`
`reg`, `total`, `item`, `times`	sums: the cash register’s digits, starting total, item and count	`4`, `1000`, `0.4`, `10`
`adds`	sums: how many times 0.1 is added (up to 2 × 10⁷)	`10000000`
`seed`	sums: which random list of 1,000 floats	`1`
`deg`	Horner: the degree n of (x − 1)ⁿ (1–24)	`10`
`at`	Horner: the x evaluated in the table	`1.02`
`w`	Horner: the plot shows 1 ± w	`0.05`
`start`	recurrence: N for the backward run from I_N = 0	`40`
`rprec`	recurrence: precision of the forward run, `f64` or `f32`	`f32`
`eps`	conditioning: ε in x + (1 + ε)y = 2 + ε	`0.001`
`nudge`	conditioning: change to the second right-hand side	`0.001`
`piv`	pivoting: log₁₀ of the pivot (−30 to 0)	`-20`
`nodes`	Runge: number of nodes, 2–25 (degree = nodes − 1)	`11`
`spacing`	Runge: `eq` or `cheb`	`cheb`
`interp`	Runge: `poly` or `spline`	`spline`
`example`	a book example (below); its settings win over any given beside it	`ch32-pivot`

fn

derivatives: sin, exp, cube (x³) or sqrt

sin

x

derivatives: the point (default 1, or 0.01 for √x)

1

prec

derivatives: f16, f32 or f64

f32

h

derivatives: log₁₀ of the step h, −20 to 0

-8

a, b, c

quadratic: coefficients of ax² + bx + c

1, -1e8, 1

digits

quadratic: significant digits, 16, 7 or 4

4

n

(1 + 1/n)ⁿ: the value of n

1e15

reg, total, item, times

sums: the cash register’s digits, starting total, item and count

4, 1000, 0.4, 10

adds

sums: how many times 0.1 is added (up to 2 × 10⁷)

10000000

seed

sums: which random list of 1,000 floats

1

deg

Horner: the degree n of (x − 1)ⁿ (1–24)

10

at

Horner: the x evaluated in the table

1.02

w

Horner: the plot shows 1 ± w

0.05

start

recurrence: N for the backward run from I_N = 0

40

rprec

recurrence: precision of the forward run, f64 or f32

f32

eps

conditioning: ε in x + (1 + ε)y = 2 + ε

0.001

nudge

conditioning: change to the second right-hand side

0.001

piv

pivoting: log₁₀ of the pivot (−30 to 0)

-20

nodes

Runge: number of nodes, 2–25 (degree = nodes − 1)

11

spacing

Runge: eq or cheb

cheb

interp

Runge: poly or spline

spline

example

a book example (below); its settings win over any given beside it

ch32-pivot

Book examples, as ?example=<name> (the tab is chosen for you):

Name Tab Shows

Name	Tab	Shows
`ch01-sum-order`	sum	Floating-point mode: 1,000 floats summed in different orders
`ch02-tenths`	sum	Add 0.1 again and again: the float drifts
`ch06-ill-conditioned`	cond	Nearly parallel lines: nudge b by 0.001
`ch07-horner`	horner	Horner versus the naive loop on expanded (x − 1)¹⁰
`ch08-cancellation`	quad	x² − 10⁸x + 1: the textbook small root is 25% wrong
`ch08-by-hand`	quad	The 4-digit calculator: x² + 62.10x + 1
`ch08-precision-lost`	quad	How much precision do you lose? x² + bx + 1
`ch16-difference-quotient`	diff	Difference quotient of sin at 1: h down to 10⁻²⁰
`ch16-fp-limits`	diff	Floating-point limits: x + h becomes x
`ch16-e-limit`	limit	(1 + 1/n)ⁿ drifts, then collapses to 1
`ch17-sweet-spot`	diff	Finding the sweet spot (64-bit)
`ch17-sweet-spot-f32`	diff	The sweet spot moves to larger h in 32-bit floats
`ch32-cancellation`	quad	x² + 10⁸x + 1 and the stable rewrite
`ch32-condition`	cond	x + 1.0001y = 2.0001: κ ≈ 4 × 10⁴
`ch32-recurrence`	recur	Iₙ = 1 − n·Iₙ₋₁: forwards explodes, backwards converges
`ch32-v-curve`	diff	The V-shaped error: forward and central differences
`ch32-halving-h`	diff	Halving h in float16, float32 and float64
`ch32-complex-step`	diff	Complex step: no subtraction, h = 10⁻²⁰ works
`ch32-pivot`	pivot	A 10⁻²⁰ pivot, with and without partial pivoting
`ch32-kahan`	sum	Kahan’s cash register: 1000 + ten × 0.4 on 4 digits
`ch32-tenths`	sum	0.1 added ten million times: naive versus Kahan
`ch32-runge`	runge	Runge’s phenomenon: degree 10, equally spaced
`ch32-runge-chebyshev`	runge	Degree 10 at Chebyshev nodes
`ch32-runge-spline`	runge	A cubic spline through 11 equally spaced points

ch01-sum-order

sum

Floating-point mode: 1,000 floats summed in different orders

ch02-tenths

sum

Add 0.1 again and again: the float drifts

ch06-ill-conditioned

cond

Nearly parallel lines: nudge b by 0.001

ch07-horner

horner

Horner versus the naive loop on expanded (x − 1)¹⁰

ch08-cancellation

quad

x² − 10⁸x + 1: the textbook small root is 25% wrong

ch08-by-hand

quad

The 4-digit calculator: x² + 62.10x + 1

ch08-precision-lost

quad

How much precision do you lose? x² + bx + 1

ch16-difference-quotient

diff

Difference quotient of sin at 1: h down to 10⁻²⁰

ch16-fp-limits

diff

Floating-point limits: x + h becomes x

ch16-e-limit

limit

(1 + 1/n)ⁿ drifts, then collapses to 1

ch17-sweet-spot

diff

Finding the sweet spot (64-bit)

ch17-sweet-spot-f32

diff

The sweet spot moves to larger h in 32-bit floats

ch32-cancellation

quad

x² + 10⁸x + 1 and the stable rewrite

ch32-condition

cond

x + 1.0001y = 2.0001: κ ≈ 4 × 10⁴

ch32-recurrence

recur

Iₙ = 1 − n·Iₙ₋₁: forwards explodes, backwards converges

ch32-v-curve

diff

The V-shaped error: forward and central differences

ch32-halving-h

diff

Halving h in float16, float32 and float64

ch32-complex-step

diff

Complex step: no subtraction, h = 10⁻²⁰ works

ch32-pivot

pivot

A 10⁻²⁰ pivot, with and without partial pivoting

ch32-kahan

sum

Kahan’s cash register: 1000 + ten × 0.4 on 4 digits

ch32-tenths

sum

0.1 added ten million times: naive versus Kahan

ch32-runge

runge

Runge’s phenomenon: degree 10, equally spaced

ch32-runge-chebyshev

runge

Degree 10 at Chebyshev nodes

ch32-runge-spline

runge

A cubic spline through 11 equally spaced points