Floating point and testing numerical code

Numerical code can be wrong in two ways. It can have a bug, or it can be a correct formula that the computer can't evaluate accurately. The second kind is harder to spot, because the code looks right and usually gives sensible answers. This lesson explains how computers store numbers, shows a textbook formula failing on realistic data, and sets out how to test numerical code so that both kinds of error get caught.

How doubles work

A double stores a number as $\pm m \times 2^{e}$, with 53 bits for the significand $m$. That gives about 16 significant decimal digits, whatever the size of the number. The consequence that matters: the gaps between representable numbers grow with the numbers.

Near 1, adjacent doubles are $2.2 \times 10^{-16}$ apart. Near $10^9$ they are about $1.2 \times 10^{-7}$ apart. Near $10^{16}$ they are 2 apart, so $10^{16} + 1$ can't be stored: it rounds back to $10^{16}$.

Most decimal fractions can't be stored exactly either. 0.1 in binary is a repeating fraction, like 1/3 in decimal, so 0.1 + 0.2 gives 0.30000000000000004. That last digit is one ulp: the result is the nearest double to the true sum of two slightly-off inputs, which isn't the double nearest 0.3.

Each individual operation is as accurate as it can be: the result is the exact answer, correctly rounded. Problems come from what happens to those tiny rounding errors over a whole calculation.

Catastrophic cancellation

The sample variance has a well-known shortcut that needs only one pass over the data:

It is algebraically identical to the two-pass definition $\frac{1}{n-1}\sum (x_i - \bar{x})^2$. Numerically it is a trap. When the values are large compared with their spread, like prices of around 10,000 that move by a few units, both sums are huge and almost equal. Subtracting them leaves only the digits where they differ, and those are exactly the digits rounding already damaged.

Below, the same 200 numbers with spread about 1 are shifted by an offset. Shifting never changes the variance. Watch what the three methods report:

At an offset of ten thousand all three agree to seven digits. By a million the shortcut is off in the fourth digit, by ten million it is off by more than 10%, and around $10^8$ it returns nonsense: a negative variance, which is impossible, or one a thousand times too big. The reason is the readout above the chart. When $\sum x_i^2$ is about $2 \times 10^{18}$, adjacent doubles are 256 apart, while the quantity being computed, the sum of squared deviations, is about 212. There is nothing left to compute it from.

The two-pass method subtracts the mean before squaring, so it works with small numbers throughout. Welford's algorithm does the same thing in a single pass, updating a running mean and a running sum of squared deviations, which also makes it right for streaming data such as live risk.

Key idea. Subtracting nearly equal numbers destroys precision. If a formula subtracts two large quantities to get a small one, find a version that works with the small differences directly.

Adding many numbers

Addition loses information too. Adding 1 to $10^{16}$ does nothing, because the result rounds back. Summing a long list of small numbers into a large running total loses a little on every step. The fixes:

• Compensated summation (Kahan, or Neumaier's variant) carries the lost low-order bits separately and adds them back. Python's math.fsum is exact, and since Python 3.12 the built-in sum compensates for floats too. NumPy's sum uses pairwise summation, which keeps the error small without being exact.
• Work in sensible units. Returns rather than price levels, or prices relative to a reference, keep numbers near 1 where doubles are densest.

In code

import math

import numpy as np

print(0.1 + 0.2 == 0.3, 0.1 + 0.2)
print(f"gap between doubles near 1: {math.ulp(1.0):.3g}, near 1e9: {math.ulp(1e9):.3g}")

rng = np.random.default_rng(3)
noise = rng.normal(0.0, 1.0, size=200)
true_var = noise.var(ddof=1)

def naive_var(x):
    n = len(x)
    return (np.sum(x * x) - np.sum(x) ** 2 / n) / (n - 1)

def two_pass_var(x):
    return np.sum((x - x.mean()) ** 2) / (len(x) - 1)

for offset in (0.0, 1e4, 1e6, 1e7, 1e8, 1e9):
    x = noise + offset
    print(f"offset {offset:>8.0e}: naive {naive_var(x):>10.4f}   two-pass {two_pass_var(x):.4f}   true {true_var:.4f}")

# Adding a small number to a huge one loses it; math.fsum keeps the lost low-order bits.
values = [1e16, 1.0, -1e16]
total = 0.0
for v in values:
    total += v
print(total, math.fsum(values))

# Testing numerical code: compare with a tolerance, and test properties.
assert not (0.1 + 0.2 == 0.3)
assert math.isclose(0.1 + 0.2, 0.3, rel_tol=1e-12)
np.testing.assert_allclose(two_pass_var(noise + 1e9), true_var, rtol=1e-6)
print("tests passed")

It prints False 0.30000000000000004 and gaps of 2.22e-16 near 1 and 1.19e-07 near $10^9$. The shortcut agrees with the true variance of 1.0664 up to an offset of $10^6$, gives 1.0653 at $10^7$, and 0.0000 from $10^8$ on, while the two-pass method stays at 1.0664 throughout. The plain loop over [1e16, 1.0, -1e16] returns 0.0; math.fsum returns the correct 1.0.

Testing numerical code

Numerical code needs tests as much as any other code, but == is the wrong tool for floats. What works:

• Compare with a tolerance. math.isclose(a, b, rel_tol=...) or np.testing.assert_allclose(actual, expected, rtol=..., atol=...). Use a relative tolerance for values far from zero and an absolute one near zero, where relative error is meaningless. Choose the tolerance from what the method can achieve, not from what makes the test pass.
• Test against known answers. A closed-form case, a hand calculation, or a value from an independent implementation. The Monte Carlo pricer is tested against Black–Scholes for exactly this reason.
• Test properties. The variance shouldn't change when you shift the data; a portfolio's weights should sum to 1; a call price should rise with volatility. Property tests catch the cancellation bug above even when no reference value exists.
• Test where things get hard. Large offsets, tiny values, empty or single-element inputs, and values that differ only in the last digits are where numerical bugs live.
• Randomised tests need fixed seeds and honest tolerances. For a simulation, work out the standard error and allow a few of them. A test that fails one run in twenty, or that can never fail, is worse than none.

Key idea. Never test floats with ==. Test against independent references with a justified tolerance, and test properties that must hold whatever the data.

Where this shows up in quant work

• Prices and P&L. Summing millions of small P&L numbers, or computing variance from raw price levels, are exactly the patterns above. Risk systems use stable, often streaming, algorithms.
• Money. Cash balances and order quantities are usually stored as integers (cents, lots, ticks) or decimals, never as binary floats, so that ledgers reconcile to the cent.
• Model validation. Validation teams compare implementations against independent ones and check properties like no-arbitrage bounds: the testing ideas above, applied systematically.

Exercises

• Find the smallest power of 10, $10^k$, for which 10.0**k + 1.0 == 10.0**k in Python. Explain the answer using the gap between doubles.
• Implement Welford's algorithm in Python and check it against np.var(x, ddof=1) on the shifted data from the code above.
• Write a property test for a function that computes portfolio volatility from weights and a covariance matrix. Which properties hold for every valid input?

Key takeaways

• Doubles carry about 16 significant digits; the gap between them grows with their size.
• Subtracting nearly equal large numbers (catastrophic cancellation) destroys precision. Use two-pass or Welford for variance.
• Long sums lose small terms; use compensated or pairwise summation, and keep numbers near 1.
• Test floats with justified tolerances, known answers and properties, never with ==.