Floating-Point Numbers and Stable Computation

Why does 0.1 + 0.2 differ from 0.3? Why can rewriting the same mathematical formula make its computed result more accurate? The answers depend on how numbers are stored, how each operation rounds, and how sensitive the problem is to its inputs. The numerical-analysis overview introduces the main kinds of error; this note develops them through actual calculations.

Sign, exponent and significand

Goldberg’s floating-point paper, reprinted by permission by Oracle, describes a floating-point number using a sign, a significand and a power of the base. For a normal IEEE 754 binary64 value,

$$x=(-1)^s(1+f)2^e,\qquad 0\le f<1.$$

Of the 64 bits, one stores the sign, 11 store the exponent and 52 store the fractional part. The leading binary digit of a normal value is implicitly 1, giving 53 bits of precision. If the exponent field holds the integer $E$, then $e=E-1023$. The fields $E=0$ and $E=2047$ are reserved for zeros, subnormals, infinities and NaNs; the normal-value formula does not apply to them.

Precision is not a fixed number of places after the decimal point. Within $[2^e,2^{e+1})$, adjacent positive normal values are separated by

$$\Delta=2^{e-52}.$$

Doubling the magnitude doubles the spacing. The upward gap at 1 is $2^{-52}$; at 2 it is $2^{-51}$. Near $10^{16}$ the gap is already 2, so adding 1 may lose the increment. Smaller floating-point formats change precision and range. The low-bit weight approximations discussed in model quantization also involve scales and encodings; total bit count alone does not determine computational precision.

Decimal display and stored values

Python’s floating-point tutorial explains why the binary expansion of 0.1 repeats forever and why Python usually displays the shortest decimal representation that recovers the same stored value. Displaying 0.1 does not mean storing exactly $1/10$. Printing more digits does not change the stored value either.

The examples use the standard library of CPython 3.9 or later (math.ulp was added in 3.9), assuming binary64 float and round-to-nearest arithmetic with ties to even. The first assertion checks the binary base and precision. Each code block can run independently.

import math
import sys
from fractions import Fraction

assert sys.float_info.radix == 2 and sys.float_info.mant_dig == 53
x = 0.1
print(repr(x))
print(format(x, ".17g"))
print(x.as_integer_ratio())
print(0.1 + 0.2, (0.1 + 0.2) == 0.3)
error = Fraction.from_float(x) - Fraction(1, 10)
print(error)
print(float(error), float(error / Fraction(1, 10)))
print(error / Fraction.from_float(math.ulp(x)))
print(math.ulp(1.0), math.ulp(2.0), math.ulp(1e16))
print(1e16 + 1.0 == 1e16)

0.1
0.10000000000000001
(3602879701896397, 36028797018963968)
0.30000000000000004 False
1/180143985094819840
5.551115123125783e-18 5.551115123125783e-17
2/5
2.220446049250313e-16 4.440892098500626e-16 2.0
True

as_integer_ratio() reveals the exact stored fraction. Subtracting $1/10$ from that fraction gives $1/180143985094819840$, about $5.55\times10^{-18}$. Arithmetic can introduce further rounding: 0.1 + 0.2 ends up one representable step above the stored 0.3.

To retain the meaning of decimal input, construct Decimal from a string; for exact rational values, use Fraction. Converting to either representation after first converting to float preserves the rounded value, not the original input.

Absolute error, relative error and ulps

For the true value $y$ and computed value $\hat y$,

$$E_{\mathrm{abs}}=|\hat y-y|,\qquad E_{\mathrm{rel}}=\frac{|\hat y-y|}{|y|}\quad(y\ne0).$$

Absolute error has the same units as the result; relative error measures the fraction lost. Relative error is undefined when the true value is zero. Near zero, a small absolute error can still be a large relative error.

An ulp, or unit in the last place, measures error using the value of the significand’s last digit. The representation error of 0.1 above is exactly $2/5$ of math.ulp(0.1). For a positive finite value, math.ulp gives the upward neighbouring gap, with a separate rule for the largest finite value. At powers of two, the downward and upward gaps can differ, so specify the reference when measuring ulps.

For normal values, correct rounding to nearest introduces at most half a local ulp of error. For a nonzero normal result of a basic operation, without overflow or underflow, the usual model is

$$\operatorname{fl}(a\circ b)=(a\circ b)(1+\delta), \qquad |\delta|\le u=2^{-53},$$

where $\circ$ denotes a basic arithmetic operation. This unit roundoff $u$ is half of sys.float_info.epsilon, the upward gap $2^{-52}$ at 1. A bound for one operation does not bound the relative error of an entire program by $u$.

Rounding, overflow, underflow and special values

Rounding to nearest chooses the representable value closest to the exact result. At a midpoint, it chooses the value whose last significand digit is even. 1e16 + 1 is such a midpoint: the earlier trace shows it rounding back to 1e16.

Situation	Meaning in binary64	Effect on computation
Overflow	The result exceeds the finite range; the largest finite value is about $1.798\times10^{308}$	With rounding to nearest, sufficiently large results become signed infinity; a language interface may raise an exception instead
Subnormal value	Satisfies $0<\lvert x\rvert<2^{-1022}$ (magnitude below the smallest positive normal value), represented with an implicit leading 0	Gradual underflow retains small values with fixed spacing $2^{-1074}$, while relative precision decreases
Rounding to zero	Too small to remain a nonzero subnormal	Half the smallest positive subnormal rounds to zero under ties to even
Infinity	inf and -inf are special values	They can represent range overflow or some unbounded results; inf - inf produces NaN
NaN	“Not a number”, which can arise from an invalid operation	Unequal even to itself; test with math.isnan, and use math.isfinite to check finiteness

A subnormal value need not be inexact. Zero has both positive and negative encodings. Goldberg’s historical paper explains gradual underflow and special values, while actual language interfaces need their own definitions. Python’s math exception conventions specify exceptions for some invalid operations and overflows rather than returning IEEE special values directly:

import math
import sys

print(sys.float_info.min)
tiny = math.ulp(0.0)
print(tiny, tiny / 2)
print(1e308 * 1e308)
nan = float("inf") - float("inf")
print(nan, nan == nan, math.isnan(nan))
for operation in (
    lambda: 1.0 / 0.0,
    lambda: math.sqrt(-1.0),
    lambda: math.exp(1000.0),
):
    try:
        operation()
    except (ZeroDivisionError, ValueError, OverflowError) as exc:
        print(type(exc).__name__)

2.2250738585072014e-308
5e-324 0.0
inf
nan False True
ZeroDivisionError
ValueError
OverflowError

5e-324 is the short decimal display of the smallest positive subnormal, not its exact decimal value. Scaling intermediate quantities into a suitable range often helps more than trying to recover after they become zero or infinite. A finite final result alone does not establish accuracy.

Cancellation: rewrite an equivalent formula

Subtracting nearby values can make their existing errors dominate the result. The subtraction itself may even be exact; information may already have been lost while computing its operands. Goldberg’s section on cancellation distinguishes this case from subtraction of exact operands.

Consider $g(x)=\sqrt{1+x}-1$ near zero. For real $x\ge-1$, multiplying by the conjugate gives

$$g(x)=\frac{(\sqrt{1+x}-1)(\sqrt{1+x}+1)}{\sqrt{1+x}+1} =\frac{x}{\sqrt{1+x}+1}.$$

With $x=10^{-16}$, the direct expression first rounds 1 + x to 1, then subtracts 1 to return zero. The rearranged expression keeps the small quantity $x$ in the numerator. Its denominator is near 2, so even rounding that denominator to 2 preserves a result of about $5\times10^{-17}$.

import math
from decimal import Decimal, localcontext

x = 1e-16
naive = math.sqrt(1.0 + x) - 1.0
stable = x / (math.sqrt(1.0 + x) + 1.0)
with localcontext() as ctx:
    ctx.prec = 80
    dx = Decimal.from_float(x)
    reference = dx / ((Decimal(1) + dx).sqrt() + Decimal(1))
    print(f"reference = {reference:.18e}")
    for name, value in [("naive", naive), ("stable", stable)]:
        absolute = abs(Decimal.from_float(value) - reference)
        if reference != 0:
            relative = absolute / abs(reference)
            print(f"{name} = {value:.18e}; relative error = {relative:.3e}")
        else:
            print(f"{name} = {value:.18e}; absolute error = {absolute:.3e}")

reference = 4.999999999999999770e-17
naive = 0.000000000000000000e+00; relative error = 1.000e+0
stable = 4.999999999999999895e-17; relative error = 2.500e-17

The reference uses 80 decimal digits and is a high-precision approximation. Decimal.from_float(x) preserves the same binary input, so this comparison measures algorithmic error without mixing in the difference between decimal and binary inputs. The direct expression has relative error 1; the rearranged expression has relative error about $2.5\times10^{-17}$. Independently, the expansion $g(x)=x/2-x^2/8+O(x^3)$ gives about $5\times10^{-17}$, with a second term of about $-1.25\times10^{-33}$.

Setting x to 0 also makes the true result zero, so the code reports absolute error instead. The rearranged expression avoids cancellation near zero, but not underflow of the result itself: for the smallest positive subnormal input, the result is about half that input and still rounds to zero in binary64.

Specialized library functions use the same principle: log1p(x) computes $\log(1+x)$ accurately near zero, and expm1(x) avoids the cancellation in directly evaluating $e^x-1$. Their real-valued domains still apply. Such rewrites also affect the intermediate values traversed by automatic differentiation; autodiff does not itself remove floating-point error.

Summation order and compensated summation

Real addition is associative; floating-point addition generally is not. Take $10^{16}$, ten copies of 1 and $-10^{16}$. The exact sum is 10. Adding sequentially from left to right can round away every 1 in the large partial sum; cancelling the two large values first retains them.

Goldberg’s summation-error analysis explains Kahan compensated summation: alongside the partial sum, retain low-order information lost in one addition to correct the next. The sequential loop below is explicit so that it does not depend on the implementation of built-in sum in a particular Python version.

import math

def sequential(values):
    total = 0.0
    for value in values:
        total += value
    return total

def kahan(values):
    total = correction = 0.0
    for value in values:
        adjusted = value - correction
        new_total = total + adjusted
        correction = (new_total - total) - adjusted
        total = new_total
    return total

values = [1e16] + [1.0] * 10 + [-1e16]
print(sequential(values))
print(sequential([1e16, -1e16] + [1.0] * 10))
print(kahan(values), math.fsum(values))
hard = [1e16, 1.0, -1e16]
print(kahan(hard), math.fsum(hard))

0.0
10.0
10.0 10.0
0.0 1.0

correction = (new_total - total) - adjusted estimates the discrepancy between this addition and the intended increment; the next step subtracts that correction. It recovers the lost contributions in the ten-ones example. But the last line shows ordinary Kahan still losing 1 for [1e16, 1.0, -1e16]. Compensation is not exact arithmetic.

Python’s math.fsum tracks multiple intermediate partial sums and returns the exact sums in both examples. Its accuracy relies on IEEE arithmetic guarantees and the usual ties-to-even rounding mode. The documentation also notes that extended-precision addition on some non-Windows builds can double-round an intermediate value and affect the last bit.

For $n\ge1$ finite floating-point inputs, treat their stored values as exact and let $S=\sum_i x_i$. The following sequential-summation bound counts addition rounding, excluding input conversion error. Assume no overflow or underflow, nonzero results satisfy the earlier single-operation model, and exactly zero additions have $\delta=0$. An input can pass through up to $n-1$ rounding steps; expanding their multiplicative factors gives the absolute error bound

$$|\hat S-S|\le\gamma_{n-1}\sum_{i=1}^n|x_i|, \qquad \gamma_k=\frac{ku}{1-ku}\quad(ku<1).$$

When $ku$ is small, $\gamma_k$ is approximately $ku$. If $S\ne0$, the relative error bound also includes the factor $\sum_i|x_i|/|S|$, which can be large when positive and negative terms nearly cancel.

For large datasets, balanced pairwise summation organizes additions as a tree, reducing the longest chain from $n-1$ to $\lceil\log_2 n\rceil$. Under the same rounding conditions, the bound replaces $\gamma_{n-1}$ with $\gamma_k$ for that tree depth. Adding small magnitudes first often reduces their chance of disappearing into a large partial sum. With mixed signs, however, no simple ordering guarantees the best result. Changing parallel reduction order can also change the last bits. Choose ordering, compensation or higher precision according to the magnitudes and required error. The NumPy note covers array reductions and comparisons.

Conditioning, backward stability and tolerances

Conditioning measures a mathematical problem’s sensitivity to input perturbations. For $y=a-b\ne0$, if each input has a relative perturbation of at most $\varepsilon$, the triangle inequality gives

$$|\Delta y|\le\varepsilon(|a|+|b|),\qquad \frac{|\Delta y|}{|y|}\le\kappa\varepsilon, \qquad \kappa=\frac{|a|+|b|}{|a-b|}.$$

Treating $a=1.000001$ and $b=1.000000$ as exact real inputs, their difference is $10^{-6}$ and $\kappa=2000001$. Relative input uncertainty of $10^{-8}$ gives an absolute output uncertainty bound of $2.000001\times10^{-8}$, or about 2% relative uncertainty. Even exact computation cannot remove that input uncertainty.

Backward stability asks whether the computed result is the exact answer for nearby inputs. If computing $f(x)$ gives $\hat y=f(x+\Delta x)$ and the input perturbation is small in a specified norm or componentwise scale, the backward error is small. A backward-stable algorithm must provide perturbation bounds commensurate with rounding precision throughout its applicable input range; finding an arbitrary nearby input for a single result is insufficient. Forward error directly compares $\hat y$ with $f(x)$. For small perturbations, differentiable $f$ and $f(x)\ne0$, the first-order relation is roughly “relative forward error $\lesssim$ condition number × relative backward error”. An ill-conditioned problem can have substantial forward error even with a backward-stable algorithm.

The square-root example has a different cause. For $x>0$, its relative condition number is

$$\left|\frac{xg'(x)}{g(x)}\right| =\frac{\sqrt{1+x}+1}{2\sqrt{1+x}} \longrightarrow1\quad(x\to0^+).$$

The problem is well-conditioned here, yet the direct expression loses the result. This distinction helps decide whether to improve the inputs or rewrite the algorithm.

Tolerances should follow an error budget: input uncertainty, discretization or approximation error, and algorithmic rounding error. An absolute tolerance needs the result’s units; a relative tolerance needs a meaningful nonzero scale. Machine precision describes representation and cannot replace that budget.

For finite values, math.isclose uses

$$|a-b|\le\max\bigl(\text{rel\_tol}\max(|a|,|b|),\text{abs\_tol}\bigr).$$

import math

print(math.isclose(0.1 + 0.2, 0.3, rel_tol=0.0, abs_tol=math.ulp(0.3)))
print(math.isclose(1e-12, 0.0, rel_tol=1e-9))
print(math.isclose(1e-12, 0.0, rel_tol=1e-9, abs_tol=1e-11))

True
False
True

The first line allows one ulp referenced to the stored 0.3, suitable for this short, analysed calculation. The next two lines show the need for a meaningful absolute tolerance near zero. 1e-11 is a demonstration threshold, not a general default. NaN is never close to anything; infinities are close only to the same signed infinity. For arrays, choose tolerances using NumPy’s interface definition rather than applying the math.isclose formula to another API.

In practice, establish input and output scales, look for nonfinite intermediates and cancellation, choose reliable summation or specialized functions, then check error against an exact solution, a high-precision reference or a proven bound. Agreement between two algorithms shows consistency; accuracy needs independent evidence.

My reading · Opened notes are remembered in this browser.