9  Floating-Point Numbers, Error, and Stability

9.1 Learning objectives

After completing this unit, you will be able to:

• explain the basic model of binary floating-point numbers;
• interpret ulp, machine epsilon, and the overflow and underflow limits;
• recognize cancellation that destroys significant digits;
• distinguish forward error from backward error;
• distinguish problem conditioning from algorithmic stability;
• choose more stable reformulations and summation methods;
• compare naive evaluation with scipy.special on a stated domain;
• set tolerances as an error budget before inspecting the results; and
• limit experimental evidence to the cases and environments actually run.

Local prerequisites: Units 1, 2, and 3; integer powers, square roots, absolute values, rational functions, and relative error at A30 level. No calculus is required.

9.2 A finite binary model

In this project’s environment, float uses a 64-bit binary format. In simplified form, a finite normal number is stored as

$$(-1)^s\,(1.b_1b_2\ldots b_{52})_2\,2^e.$$

There is one sign bit, a bounded exponent, and 53 bits of significand precision when the leading one before the binary point is counted. Because there are only finitely many bit patterns, only finitely many real numbers can be stored. Operation results are usually rounded to the nearest representable number.

The hex() method reveals the binary value of a float without losing information:

x = 0.1
(x, x.hex(), x.as_integer_ratio())

Listing 9.1

(0.1, '0x1.999999999999ap-4', (3602879701896397, 36028797018963968))

The fraction returned by as_integer_ratio() is the rational value actually stored. It is close to $1/10$, but it is not $1/10$.

9.3 ulp, epsilon, and changing spacing

One unit in the last place (ulp) is the local spacing between adjacent representable numbers. This spacing changes with the magnitude of the number.

import math
import sys

(
    sys.float_info.epsilon,
    math.ulp(1.0),
    math.ulp(float(2**53)),
)

Listing 9.2

(2.220446049250313e-16, 2.220446049250313e-16, 2.0)

Around 1, ulp(1.0) equals sys.float_info.epsilon: the difference between 1 and the next larger float. Some books use “machine epsilon” for this number. Others use the unit roundoff $u=\varepsilon/2$ in the round-to-nearest model. Because terminology varies, state the definition you use.

Around $2^{53}$, one ulp is 2. Consequently:

large_value = float(2**53)
(large_value + 1.0 == large_value, large_value + 2.0 == large_value + 2)

Listing 9.3

(True, True)

This is not a mistake in Python’s integer addition. The issue is that the result is being forced back onto the float64 grid, whose spacing at this scale is 2.

9.4 Overflow, subnormals, and underflow

Precision and range are two different limits. sys.float_info.max is the largest finite float. Multiplying it by 2 produces positive infinity in the tested environment. At the small end, sys.float_info.min is the smallest positive normal number, not the smallest positive number.

Subnormal numbers fill part of the gap toward zero with reduced precision. math.ulp(0.0) gives the smallest positive subnormal. Dividing it by two gives zero because there is no smaller positive representable number.

largest_float = sys.float_info.max
smallest_subnormal = math.ulp(0.0)
(
    largest_float * 2.0,
    sys.float_info.min,
    smallest_subnormal,
    smallest_subnormal / 2.0,
)

Listing 9.4

(inf, 2.2250738585072014e-308, 5e-324, 0.0)

Check values with math.isfinite before passing an infinite result to the next stage, and check separately for zero caused by underflow. Checks after the event do not replace an analysis of scale before computation.

9.5 Cancellation and stable reformulation

Subtracting two nearly equal numbers can eliminate most of their significant digits. For large $x$, consider

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

At $x=10^{16}$, the addition x + 1.0 already rounds back to x. Direct subtraction then produces zero. Before computing, rationalize algebraically:

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

x = 1e16
direct = math.sqrt(x + 1.0) - math.sqrt(x)
stable = 1.0 / (math.sqrt(x + 1.0) + math.sqrt(x))
(x + 1.0 == x, direct, stable)

Listing 9.5

(True, 0.0, 5e-09)

The second form avoids subtracting two nearly equal square roots. This reformulation follows from an algebraic identity over the real numbers; the Unit 6 program compares the two implementations with a high-precision Decimal reference. A high-precision reference remains a checking tool, not a proof in its own right.

9.6 SciPy: stable exprel near zero

SciPy provides a specialized implementation of the function

$$E(x)=\frac{e^x-1}{x}\quad (x\ne0), \qquad E(0)=1.$$

Its mathematical domain is all real numbers once the value at zero is defined by this continuous extension. This unit’s detailed lab uses the narrower computational domain $0<|x|\le10^{-4}$ because the reported relative backward error divides by $|x|$. The Decimal oracle does not compute exp(x) - 1: it sums the power series for $E$ and $E'$ without cancellation. Its precision is at least 80 digits and increases with the exponent of $x$, retaining 50 guard digits beyond the order of $x$. The first correction $x/2$ therefore remains resolved even for the smallest binary64 subnormal. The upper bound and the exclusion of $x=0$ delimit the lab’s error report, not the mathematical domain of exprel.

Direct evaluation computes exp(x), subtracts 1, and divides by x. When x is small, exp(x) and 1 are nearly equal, so subtraction can discard significant digits. For more extreme values, exp(x) may round to exactly 1, making the naive numerator zero.

import math
import scipy
from scipy import special

x = 1e-8
naive = (math.exp(x) - 1.0) / x
stable = float(special.exprel(x))

extreme_x = 1e-16
extreme_naive = (math.exp(extreme_x) - 1.0) / extreme_x
extreme_stable = float(special.exprel(extreme_x))

(scipy.__version__, naive, stable, extreme_naive, extreme_stable)

Listing 9.6

('1.15.2', 0.999999993922529, 1.000000005, 0.0, 1.0)

This difference does not automatically mean that the mathematical problem is ill-conditioned. A measure of local relative conditioning is

$$\kappa_E(x)=\left|\frac{xE'(x)}{E(x)}\right|.$$

For $x\ne0$,

$$E'(x)=\frac{(x-1)e^x+1}{x^2}.$$

Near zero, $E(x)$ approaches 1, $E'(x)$ approaches $1/2$, and $\kappa_E(x)$ is approximately $|x|/2$. This derivative formula is supplied as a diagnostic tool; deriving it is not an A30 prerequisite. At $x=10^{-8}$, the condition number is about $5\times10^{-9}$, so small relative perturbations in the input should not be amplified. The large error in naive evaluation is therefore algorithmic instability, not inherent sensitivity of the function.

The companion script compares both results with an adaptive-precision Decimal series reference on the stated lab domain. The reported forward error is

$$e_f=|\widehat E-E(x)|.$$

To interpret the result in terms of backward error, the script uses the first-order estimate

$$|\Delta x|\approx\frac{e_f}{|E'(x)|}.$$

This asks approximately how much the input would have to change for the computed output to be the exact output for that changed input. This estimate is not an exact inverse solution, and the JSON labels it accordingly. At $x=10^{-8}$, the naive evaluation’s estimated relative backward error is greater than 1, whereas SciPy’s estimate is far below $10^{-6}$. SciPy’s relative forward error is also below $10^{-15}$ in the tested run.

These results provide evidence for the recorded values, SciPy version, and environment. They give a concrete example of a specialized implementation being more accurate than the naive formula near zero. They do not prove that scipy.special.exprel is stable for every real number or every backend. A general claim still requires analysis of the algorithm, the rounding model, and range handling.

9.7 Forward error and backward error

Suppose a problem asks for $y=f(x)$ and a program returns $\widehat y$.

• Forward error measures the distance from the desired answer: $|\widehat y-f(x)|$.
• Backward error asks how far the input must change for the output to be an exact answer: find $\Delta x$ such that $\widehat y=f(x+\Delta x)$.

For the problem $y=\sqrt a$, squaring the output gives an input for which that output would be exact. A readily checkable measure of relative backward error is therefore

$$\frac{|\widehat y^2-a|}{|a|}.$$

y_hat = math.sqrt(2.0)
residual = y_hat * y_hat - 2.0
(y_hat, residual)

Listing 9.7

(1.4142135623730951, 4.440892098500626e-16)

A small backward error means that the result is an exact solution to a slightly perturbed problem. Whether its forward error is also small depends on the conditioning of the problem.

9.8 Conditioning is not stability

Conditioning is a property of the mathematical question: how much does the output change when the input changes slightly? Algorithmic stability is a property of the method of computation: does internal rounding introduce errors far beyond the problem’s sensitivity?

For

$$f(x)=\frac{1}{1-x},$$

inputs close to 1 are highly sensitive. At $x=0.99999999$, a perturbation of about $10^{-12}$ can be amplified by roughly $10^8$ in the relative output error. This problem is ill-conditioned near the singular point $x=1$. Changing programming languages cannot remove that mathematical sensitivity.

By contrast, the problem $d(x)=\sqrt{x+1}-\sqrt{x}$ for large $x$ does not require a zero answer. Direct subtraction loses information, whereas the rationalized form retains it. Here, the main difference is the stability of the evaluation method.

Four possibilities must be distinguished:

Problem	Algorithm	Practical meaning
well-conditioned	stable	small errors are generally expected
well-conditioned	unstable	the algorithm wastes information
ill-conditioned	stable	results can still be sensitive to the data
ill-conditioned	unstable	problem sensitivity and algorithmic error compound

9.9 Summation and order of operations

Floating-point addition is not associative. Over the real numbers,

$$10^{16}+1-10^{16}=1.$$

However, a naive algorithm that sums from left to right loses the 1 before the final subtraction.

values = [1e16, 1.0, -1e16]

def naive_sum(data):
    total = 0.0
    for value in data:
        total += value
    return total

(
    naive_sum(values),
    naive_sum([1e16, -1e16, 1.0]),
    math.fsum(values),
)

Listing 9.8

(0.0, 1.0, 1.0)

math.fsum tracks the small lost parts more carefully and returns 1 for this example. The built-in sum may use an improved strategy in some Python versions; the script therefore records its result too, but does not use it as the definition of the naive algorithm. math.fsum is generally more accurate for data of mixed magnitudes, but this is not a guarantee that every summation problem is error-free. Record the method, version, data order, and magnitudes.

9.10 Tolerances as an error budget

The closeness test used in this unit has the form

$$|\widehat y-y_{ref}|\le A+R|y_{ref}|,$$

where $A$ is the absolute budget and $R$ is the relative budget. The absolute part matters near zero; the relative part scales with the reference.

The budget must come from the problem before the result is inspected. Its sources may include measurement precision, model approximation, discretization, iteration stopping criteria, and rounding. Adding the bounds for individual components gives a conservative bound when no sharper analysis is available.

isclose or a similar function only checks the rule it is given. It does not choose a tolerance, repair an algorithm, or prove that the reference is correct. Save the values of $A$ and $R$, the observed error, and the reasons for choosing them.

9.11 Deterministic records and limits of evidence

Run from the project root:

python source/code/unit06_floating.py --output output/unit06-results.json

The script records the Python and SciPy versions, the binary64 profile, cancellation, the scipy.special.exprel comparison, forward and backward errors, conditioning, summation, and error-budget decisions. Numbers that must be distinguished bit by bit are stored using repr, hex, or Decimal text. The JSON uses UTF-8, sorted keys, fixed indentation, and LF line endings. The core_sha256 field binds the record before the hash is added.

The experiment establishes the behavior of the cases actually run in the recorded environment. It can supply a counterexample to a claim such as “summation order never matters.” It does not prove a theorem about all binary64 operations.

General statements require a model, for example: each normal operation is rounded to the nearest representable value with relative error at most approximately $u$, provided no overflow or underflow occurs. Error bounds can then be derived step by step from that model. Tests check whether the implementation and environment follow the recorded assumptions.

9.12 Exercises

9.12.1 Exercise 1 - ulp and the lost increment

Around $2^{53}$, ulp is 2. Explain why float(2**53) + 1.0 == float(2**53), but adding 2 can produce the next number.

NoteHint

Imagine a grid of representable numbers whose points are spaced 2 apart.

TipAnswer and solution

At this scale, adjacent finite float values are $2^{53}, 2^{53}+2, 2^{53}+4,\ldots$. The real value $2^{53}+1$ lies exactly between two grid points. The ties-to-even rule rounds it to the value with an even significand, namely $2^{53}$. The value $2^{53}+2$ is itself representable, so the increment of 2 is not lost.

9.12.2 Exercise 2 - eliminating cancellation

Derive a stable form of $\sqrt{x+4}-\sqrt{x}$ and explain why it is better for large $x$.

NoteHint

Multiply the numerator and denominator by the sum of the two square roots.

TipAnswer and solution

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

The direct form subtracts two large, nearly equal numbers, so significant digits can be lost. The rationalized form divides 4 by a positive sum, avoiding that subtraction. The real domain considered here is $x\ge0$.

9.12.3 Exercise 3 - forward and backward error

For the approximation $\widehat y=1.414$ to $\sqrt2$, estimate the absolute forward error and the relative backward error based on the equation $y^2=2$.

NoteHint

Use $\sqrt2\approx1.41421356$ and compute $1.414^2$.

TipAnswer and solution

The absolute forward error is approximately $|1.414-1.41421356|=0.00021356$. Since $1.414^2=1.999396$, the relative backward error with respect to the input 2 is

$$\frac{|1.999396-2|}{2}=0.000302.$$

The first number compares the output with the desired square root; the second asks what relative change in the input would make that output an exact square root.

9.12.4 Exercise 4 - conditioning or algorithm?

Classify these two observations: (a) the value of $1/(1-x)$ changes greatly when $x$ near 1 is perturbed slightly; (b) subtracting two nearly equal square roots gives zero, but the rationalized form gives an accurate nonzero value.

NoteHint

Ask whether the sensitivity already belongs to the mathematical function or arises from the sequence of operations.

TipAnswer and solution

Observation (a) is primarily ill-conditioning: the function has a singular point at 1 and is inherently sensitive nearby. Even a stable algorithm cannot remove that sensitivity to the data. Observation (b) is primarily instability of direct evaluation: the same problem can be computed much more accurately after algebraic reformulation. In practice, both can occur together.

9.12.5 Exercise 5 - summation, tolerance, and limits on claims

For the data [1e16, 1.0, -1e16], the mathematical target is 1. Compare the results of naive_sum above and math.fsum. With an absolute budget of 0.5 and a relative budget of $10^{-12}$, determine which passes. Does one success prove that math.fsum is always exact?

NoteHint

Use the bound $A+R|y_{ref}|$ with $y_{ref}=1$.

TipAnswer and solution

naive_sum returns 0 because the 1 is lost when added to $10^{16}$. math.fsum returns 1 in this case. The allowed bound is $0.5+10^{-12}$; the naive algorithm’s error is 1, so it fails, whereas the error of math.fsum is 0, so it passes. This success establishes only the result of this case in the tested environment; it does not prove that math.fsum is always exact for every list.

9.12.6 SciPy exercise - stability of exprel

Run scipy_exprel_report(1e-8) from the companion code. Record the mathematical and lab domains, the relative condition number, both implementations’ relative forward errors, and their estimated relative backward errors. Determine whether the poor naive result primarily indicates an ill-conditioned problem or an unstable algorithm. Repeat the value comparison at x=1e-16 and limit your claims to the evidence actually obtained.

NoteHint

If the condition number is far below 1 but the naive forward error is large, amplification of input perturbations is not the problem. Also check whether math.exp(1e-16) == 1.0.

ImportantExecutable check

From the project root, run:

python tests/test_unit06.py Unit06FloatingTests.test_scipy_special_exprel_is_stable_near_zero -v

The check covers the lab domain, SciPy version and function, forward error, estimated backward error, conditioning, extreme cancellation, and limits of the evidence. A separate regression test also runs the smallest binary64 subnormal and checks that the oracle still captures the first correction $x/2$, rather than silently rounding it to zero.

TipAnswer and solution

Listing 9.9

import runpy

unit06 = runpy.run_path("source/code/unit06_floating.py")
report = unit06["scipy_exprel_report"](1e-8)
print("domain:", report["mathematical_domain"])
print("lab:", report["laboratory_domain"])
print("condition:", report["conditioning"]["relative_condition_number"])
print("naive forward:", report["forward_error"]["naive_relative"])
print("SciPy forward:", report["forward_error"]["scipy_relative"])
print("naive backward:", report["backward_error"]["naive_relative_estimate"])
print("SciPy backward:", report["backward_error"]["scipy_relative_estimate"])

domain: all real x
lab: 0 < |x| <= 1e-4 for binary64 inputs
condition: 5.00000000833333343795E-9
naive forward: 1.10774709322014616851E-8
SciPy forward: 4.70540214322867789081E-17
naive backward: 2.21549418274780198609E+0
SciPy backward: 9.41080427077268180010E-9

The mathematical domain is all real numbers with $E(0)=1$; the lab is restricted to $0<|x|\le10^{-4}$ for binary64 inputs. At $x=10^{-8}$, the condition number is about $5\times10^{-9}$, but the naive relative forward error exceeds $10^{-9}$ and its estimated relative backward error exceeds 1. By contrast, SciPy’s relative forward error is below $10^{-15}$ and its estimated relative backward error is below $10^{-6}$. Since the problem is well-conditioned near zero, this difference reveals cancellation in the naive algorithm.

At $x=10^{-16}$, math.exp(x) rounds to 1, so the naive formula gives 0, whereas scipy.special.exprel gives 1 in this run. These tests do not prove behavior for all inputs or backends; they check only the recorded cases and environment.

9.13 Summary

• float values form a finite binary grid; ulp changes with scale.
• Epsilon, unit roundoff, the normal range, and the subnormal range answer different questions.
• Cancellation can be avoided through algebraic reformulation before evaluation.
• scipy.special.exprel provides an executable example of a specialized function that avoids naive cancellation near zero on the stated lab domain.
• Forward error measures the answer; backward error measures the change in the problem that makes that answer exact.
• Conditioning belongs to the problem, whereas stability belongs to the algorithm.
• Summation methods and operation order affect error.
• A tolerance is a budget justified in advance, not a number adjusted to make a test pass.
• Experiments provide bounded evidence; general claims require a rounding model and a mathematical argument.