Floating-point landmarks

Floating-point behaviour is systematic: the format, the scale of the values, and the order of operations constrain what a computation can retain. This tutorial makes those constraints observable and uses them to explain the variance discrepancy from Module 1.

This activity accompanies Module 2: Understanding Floating-Point Arithmetic.

Download the executable Jupyter notebook Open in Google Colab

Learning goals

After this experiment, you should be able to:

  • inspect the floating-point format used by a Python environment;
  • compare the precision-and-range trade-offs of binary16, bfloat16, binary32, and binary64;
  • observe how representable spacing changes with magnitude;
  • predict when an increment will be rounded away;
  • classify normal, subnormal, infinite, zero, and NaN values;
  • demonstrate that evaluation order can change a computed result;
  • use local spacing to explain the Module 1 variance discrepancy.

Prerequisites

Complete the opening experiment or read Module 1: When Correct Code Produces Wrong Answers first. You need basic Python expressions and loops, but no knowledge of the IEEE 754 encoding.

Outline

  1. Confirm the floating-point format.
  2. Compare common binary formats with two revealing values.
  3. Inspect a decimal value that is not exactly representable.
  4. Map spacing at several magnitudes.
  5. Observe when an increment has no effect.
  6. Classify values at the edges of the format.
  7. Compare evaluation orders.
  8. Explain the variance discrepancy.

Confirm the working format

The published course environment uses IEEE 754 binary64 for Python’s float. The assertions make that assumption explicit: if this notebook runs in a different environment, it should stop instead of silently attaching binary64 claims to another format. Radix is another name for the base of the number system, so a radix of 2 means that this activity is using binary floating point.

import math
import struct
import sys


assert sys.float_info.radix == 2, "This activity expects binary floating point."
assert sys.float_info.mant_dig == 53, "This activity expects binary64 precision."
assert sys.float_info.max_exp == 1024, "This activity expects binary64 range."

print(f"radix:                 {sys.float_info.radix}")
print(f"significand precision: {sys.float_info.mant_dig} bits")
print(f"largest finite value:  {sys.float_info.max:.6e}")
print(f"smallest normal value: {sys.float_info.min:.6e}")
print(f"epsilon near one:      {sys.float_info.epsilon:.6e}")
radix:                 2
significand precision: 53 bits
largest finite value:  1.797693e+308
smallest normal value: 2.225074e-308
epsilon near one:      2.220446e-16

The significand precision describes how much binary detail a normal value can retain. The exponent range allows magnitudes from approximately \(10^{-308}\) to \(10^{308}\), but it does not make every number in that interval representable.

Compare common binary formats

Binary16 and bfloat16 both use 16 storage bits but allocate them differently. Binary16 has 11 bits of significand precision and 5 exponent bits; bfloat16 has 8 bits of significand precision and 8 exponent bits. Binary32 and binary64 have 24 and 53 bits of significand precision, respectively.

The helpers below emulate a declared conversion path without requiring a third-party numerical library. Python’s struct module performs binary16 and binary32 conversions. The bfloat16 helper first converts the finite input to binary32 and then rounds that bit pattern to bfloat16 using round-to-nearest, ties-to-even. It is a format-conversion experiment, not an emulator for a particular CPU, GPU, compiler, or arithmetic kernel.

def round_binary16(value):
    """Round a finite value to binary16, then return it as binary64."""
    return struct.unpack(">e", struct.pack(">e", value))[0]


def round_binary32(value):
    """Round a finite value to binary32, then return it as binary64."""
    return struct.unpack(">f", struct.pack(">f", value))[0]


def round_bfloat16_from_binary32(value):
    """Round a finite binary32 conversion to bfloat16, returning binary64."""
    if not math.isfinite(value):
        raise ValueError("This teaching helper expects a finite value.")
    binary32_bits = struct.unpack(">I", struct.pack(">f", value))[0]
    rounding_bias = 0x7FFF + ((binary32_bits >> 16) & 1)
    bfloat16_bits = (binary32_bits + rounding_bias) & 0xFFFF0000
    return struct.unpack(">f", struct.pack(">I", bfloat16_bits))[0]


format_rounders = {
    "binary16": round_binary16,
    "bfloat16": round_bfloat16_from_binary32,
    "binary32": round_binary32,
    "binary64": float,
}

Before running the next cell, make two predictions:

  1. Which formats distinguish 1.001 from 1.0?
  2. Which formats retain \(10^{-20}\) as a nonzero value?

These probes separate local precision near one from exponent range.

format_probes = {
    "1.001": 1.001,
    "1e-20": 1.0e-20,
}
format_probe_records = []

print(f"{'input':>8s}  {'format':>9s}  {'stored value':>24s}")
for input_label, input_value in format_probes.items():
    for format_name, round_value in format_rounders.items():
        stored_value = round_value(input_value)
        format_probe_records.append(
            {
                "input": input_label,
                "format": format_name,
                "stored_value": stored_value,
            }
        )
        print(f"{input_label:>8s}  {format_name:>9s}  {stored_value:24.17g}")
   input     format              stored value
   1.001   binary16              1.0009765625
   1.001   bfloat16                         1
   1.001   binary32        1.0010000467300415
   1.001   binary64        1.0009999999999999
   1e-20   binary16                         0
   1e-20   bfloat16    1.0005576689441423e-20
   1e-20   binary32    9.9999996826552254e-21
   1e-20   binary64    9.9999999999999995e-21

Binary16 changes 1.001 to its nearby value 1.0009765625, while bfloat16 rounds it to 1. Bfloat16 has coarser spacing near one because it retains fewer significand bits. Conversely, binary16 rounds \(10^{-20}\) to zero, while bfloat16 retains a normal value because its wider exponent field provides approximately the normal range of binary32.

The result is a trade-off, not a quality ranking. Whether either approximation is adequate depends on the required accuracy and on the complete computation.

Predict the first lost unit increment

If a binary format has \(p\) bits of significand precision, all consecutive integers through \(2^p\) are representable. At \(2^p\), the upward spacing is two, so the halfway sum \(2^p+1\) can round back to \(2^p\). Some larger integers remain representable, but not every integer does. Predict the threshold for each format before running the check.

significand_bits = {
    "binary16": 11,
    "bfloat16": 8,
    "binary32": 24,
    "binary64": 53,
}
integer_threshold_records = []

print(f"{'format':>9s}  {'2**p':>18s}  {'2**p + 1 changes?':>18s}")
for format_name, precision_bits in significand_bits.items():
    threshold = float(2**precision_bits)
    rounded_sum = format_rounders[format_name](threshold + 1.0)
    changed = rounded_sum != threshold
    integer_threshold_records.append(
        {
            "format": format_name,
            "significand_bits": precision_bits,
            "threshold": threshold,
            "adding_one_changes_value": changed,
        }
    )
    print(f"{format_name:>9s}  {threshold:18.0f}  {str(changed):>18s}")
   format                2**p   2**p + 1 changes?
 binary16                2048               False
 bfloat16                 256               False
 binary32            16777216               False
 binary64    9007199254740992               False

The same finite-set model predicts every row; only the format parameters change. This is why a result record should name the format instead of saying only “floating point,” “single precision,” or “16 bit.”

Inspect a stored decimal fraction

Before running the cell, predict whether the exact stored value of 0.1 is the rational number \(1/10\). The hexadecimal form is a compact, exact description of the stored binary floating-point value.

value = 0.1
numerator, denominator = value.as_integer_ratio()

print(f"usual display:       {value}")
print(f"17 significant digits: {value:.17g}")
print(f"exact stored ratio:  {numerator} / {denominator}")
print(f"exact hexadecimal:   {value.hex()}")
print(f"stored value equals 1/10: {numerator * 10 == denominator}")
usual display:       0.1
17 significant digits: 0.10000000000000001
exact stored ratio:  3602879701896397 / 36028797018963968
exact hexadecimal:   0x1.999999999999ap-4
stored value equals 1/10: False

The stored ratio is close to, but not equal to, \(1/10\). Python’s usual display uses the shortest decimal text that converts back to the same stored value. Showing more digits reveals the approximation; it does not alter it.

Map the local spacing

math.nextafter(x, math.inf) gives the next representable value above x. Predict whether the absolute gap will remain constant as x grows.

scales = [1.0e-6, 1.0, 1.0e6, float(2**53), 1.0e16]

print(f"{'x':>24s}  {'next value':>24s}  {'upward gap':>14s}")
for x in scales:
    next_value = math.nextafter(x, math.inf)
    upward_gap = next_value - x
    print(f"{x:24.17g}  {next_value:24.17g}  {upward_gap:14.6g}")
                       x                next value      upward gap
  9.9999999999999995e-07    1.0000000000000002e-06     2.11758e-22
                       1        1.0000000000000002     2.22045e-16
                 1000000        1000000.0000000001     1.16415e-10
        9007199254740992          9007199254740994               2
       10000000000000000         10000000000000002               2

The absolute gap grows with magnitude. Near one, it is approximately \(2.22\times10^{-16}\); at \(2^{53}\) and \(10^{16}\), it is 2. Values that differ by less than a local gap may map to the same floating-point value.

math.ulp(x) reports a related local spacing measure. An ulp is meaningful at the scale of x; it is not one global comparison tolerance.

Observe a lost increment

At \(2^{53}\), adjacent binary64 values are two units apart. Predict the result of adding each increment before running the cell.

base = float(2**53)
increments = [0.5, 1.0, 2.0, 3.0, 4.0]

print(f"base = {base:.0f}, upward gap = {math.ulp(base):.0f}")
print(f"{'increment':>10s}  {'computed sum':>18s}  {'changed?':>8s}")
for increment in increments:
    computed_sum = base + increment
    print(
        f"{increment:10.1f}  {computed_sum:18.0f}  "
        f"{str(computed_sum != base):>8s}"
    )
base = 9007199254740992, upward gap = 2
 increment        computed sum  changed?
       0.5    9007199254740992     False
       1.0    9007199254740992     False
       2.0    9007199254740994      True
       3.0    9007199254740996      True
       4.0    9007199254740996      True

base + 1.0 is exactly halfway between two representable values. Under the default round-to-nearest, ties-to-even rule, it rounds back to base. The increment is not zero; its contribution is absent from this stored sum.

Exercise: find an effective increment

Change test_increment and find the smallest positive integer increment for which the computed sum differs from base. Predict the answer from the local spacing before rerunning the cell.

test_increment = 1.0  # Change this value.
test_sum = base + test_increment

print(f"base:      {base:.0f}")
print(f"increment: {test_increment:.1f}")
print(f"sum:       {test_sum:.0f}")
print(f"changed:   {test_sum != base}")
base:      9007199254740992
increment: 1.0
sum:       9007199254740992
changed:   False

Check the exercise

The next cell checks the first four positive integer increments. Explain why the first effective increment agrees with the spacing observed above.

for integer_increment in range(1, 5):
    computed_sum = base + float(integer_increment)
    print(
        f"increment {integer_increment}: "
        f"sum = {computed_sum:.0f}, changed = {computed_sum != base}"
    )
increment 1: sum = 9007199254740992, changed = False
increment 2: sum = 9007199254740994, changed = True
increment 3: sum = 9007199254740996, changed = True
increment 4: sum = 9007199254740996, changed = True

Classify exceptional numerical states

The following classifier makes each state explicit. Before running it, predict which operations create infinity, a subnormal number, and NaN.

def classify_float(value):
    """Classify a Python float under the expected binary64 format."""
    if math.isnan(value):
        return "NaN"
    if math.isinf(value):
        return "+infinity" if value > 0.0 else "-infinity"
    if value == 0.0:
        sign = math.copysign(1.0, value)
        return "+zero" if sign > 0.0 else "-zero"
    if abs(value) < sys.float_info.min:
        return "subnormal"
    return "normal finite"


largest_finite = sys.float_info.max
smallest_normal = sys.float_info.min
smallest_subnormal = math.ulp(0.0)

states = [
    ("ordinary", 1.0),
    ("overflowed", largest_finite * 2.0),
    ("smallest normal", smallest_normal),
    ("half smallest normal", smallest_normal / 2.0),
    ("smallest subnormal", smallest_subnormal),
    ("underflowed", smallest_subnormal / 2.0),
    ("negative zero", -0.0),
    ("indeterminate", math.inf - math.inf),
]

print(f"{'source':>22s}  {'stored value':>24s}  {'classification':>14s}")
for source, state_value in states:
    print(
        f"{source:>22s}  {state_value:24.17g}  "
        f"{classify_float(state_value):>14s}"
    )
                source              stored value  classification
              ordinary                         1   normal finite
            overflowed                       inf       +infinity
       smallest normal   2.2250738585072014e-308   normal finite
  half smallest normal   1.1125369292536007e-308       subnormal
    smallest subnormal   4.9406564584124654e-324       subnormal
           underflowed                         0           +zero
         negative zero                        -0           -zero
         indeterminate                       nan             NaN

Subnormal values fill part of the gap between the smallest normal magnitude and zero, but with progressively fewer significant bits. The transition through that region is gradual underflow. A still smaller result rounds to signed zero.

Do not test for NaN with equality. NaN is unordered and is not equal to itself; use a classification operation such as math.isnan.

not_a_number = float("nan")

print(f"NaN equals itself:       {not_a_number == not_a_number}")
print(f"math.isnan(NaN):         {math.isnan(not_a_number)}")
print(f"positive and negative zero compare equal: {0.0 == -0.0}")
print(f"sign of +0.0:            {math.copysign(1.0, 0.0):+.0f}")
print(f"sign of -0.0:            {math.copysign(1.0, -0.0):+.0f}")
NaN equals itself:       False
math.isnan(NaN):         True
positive and negative zero compare equal: True
sign of +0.0:            +1
sign of -0.0:            -1

Change the evaluation order

Both expressions below equal 1 in real arithmetic. Predict their floating-point results by considering the first addition in each expression.

a = 1.0e16
b = -1.0e16
c = 1.0

left_grouped = (a + b) + c
right_grouped = a + (b + c)

print(f"(a + b) + c = {left_grouped:.1f}")
print(f"a + (b + c) = {right_grouped:.1f}")
print(f"results equal: {left_grouped == right_grouped}")
(a + b) + c = 1.0
a + (b + c) = 0.0
results equal: False

The first grouping cancels a and b exactly before adding c. In the second grouping, adding c to b is rounded back to b, so the later sum is zero. Floating-point addition is therefore not generally associative.

Return to the variance experiment

The shortcut formula subtracts two intermediate values near \(10^{16}\ \mathrm{ns}^2\). Inspect their local spacing and compare the result with the centered calculation and the exact reference established in Module 1.

measurements_ns = [
    100_000_004.0,
    100_000_007.0,
    100_000_013.0,
    100_000_016.0,
]

mean_ns = sum(measurements_ns) / len(measurements_ns)
centered_variance_ns2 = sum(
    (value - mean_ns) ** 2 for value in measurements_ns
) / len(measurements_ns)
mean_of_squares_ns2 = sum(
    value**2 for value in measurements_ns
) / len(measurements_ns)
square_of_mean_ns2 = mean_ns**2
shortcut_variance_ns2 = mean_of_squares_ns2 - square_of_mean_ns2

print(f"mean:                         {mean_ns:.1f} ns")
print(f"centered values:              {[value - mean_ns for value in measurements_ns]} ns")
print(f"mean of squares:              {mean_of_squares_ns2:.17g} ns²")
print(f"square of mean:               {square_of_mean_ns2:.17g} ns²")
print(f"spacing at mean of squares:   {math.ulp(mean_of_squares_ns2):.1f} ns²")
print(f"centered variance:            {centered_variance_ns2:.1f} ns²")
print(f"shortcut variance:            {shortcut_variance_ns2:.1f} ns²")
print(f"exact reference from Module 1: {22.5:.1f} ns²")
mean:                         100000010.0 ns
centered values:              [-6.0, -3.0, 3.0, 6.0] ns
mean of squares:              10000002000000122 ns²
square of mean:               10000002000000100 ns²
spacing at mean of squares:   2.0 ns²
centered variance:            22.5 ns²
shortcut variance:            22.0 ns²
exact reference from Module 1: 22.5 ns²

The centered calculation works with values near the spread of the data. The shortcut calculation first rounds much larger intermediate values whose local spacing is \(2\ \mathrm{ns}^2\), then subtracts them. The rounded intermediates no longer contain enough information to recover \(22.5\ \mathrm{ns}^2\).

This explains the observation; it does not prove that every centered variance calculation is reliable or determine an acceptable error for another use case.

Evidence record

Record the environment assumption, observations, explanation, and limitation together. This makes clear which claims may need to be revisited in another language or floating-point format.

evidence = {
    "primary_arithmetic_format": "IEEE 754 binary64",
    "rounding_assumption": "round to nearest, ties to even",
    "format_probe_records": format_probe_records,
    "integer_threshold_records": integer_threshold_records,
    "spacing_at_2**53": math.ulp(float(2**53)),
    "lost_increment": float(2**53) + 1.0 == float(2**53),
    "addition_is_associative_for_test": left_grouped == right_grouped,
    "variance_reference_ns2": 22.5,
    "centered_variance_ns2": centered_variance_ns2,
    "shortcut_variance_ns2": shortcut_variance_ns2,
    "limitation": (
        "No acceptable-error criterion or measurement uncertainty was assessed."
    ),
}
evidence
{'primary_arithmetic_format': 'IEEE 754 binary64',
 'rounding_assumption': 'round to nearest, ties to even',
 'format_probe_records': [{'input': '1.001',
   'format': 'binary16',
   'stored_value': 1.0009765625},
  {'input': '1.001', 'format': 'bfloat16', 'stored_value': 1.0},
  {'input': '1.001', 'format': 'binary32', 'stored_value': 1.0010000467300415},
  {'input': '1.001', 'format': 'binary64', 'stored_value': 1.001},
  {'input': '1e-20', 'format': 'binary16', 'stored_value': 0.0},
  {'input': '1e-20',
   'format': 'bfloat16',
   'stored_value': 1.0005576689441423e-20},
  {'input': '1e-20',
   'format': 'binary32',
   'stored_value': 9.999999682655225e-21},
  {'input': '1e-20', 'format': 'binary64', 'stored_value': 1e-20}],
 'integer_threshold_records': [{'format': 'binary16',
   'significand_bits': 11,
   'threshold': 2048.0,
   'adding_one_changes_value': False},
  {'format': 'bfloat16',
   'significand_bits': 8,
   'threshold': 256.0,
   'adding_one_changes_value': False},
  {'format': 'binary32',
   'significand_bits': 24,
   'threshold': 16777216.0,
   'adding_one_changes_value': False},
  {'format': 'binary64',
   'significand_bits': 53,
   'threshold': 9007199254740992.0,
   'adding_one_changes_value': False}],
 'spacing_at_2**53': 2.0,
 'lost_increment': True,
 'addition_is_associative_for_test': False,
 'variance_reference_ns2': 22.5,
 'centered_variance_ns2': 22.5,
 'shortcut_variance_ns2': 22.0,
 'limitation': 'No acceptable-error criterion or measurement uncertainty was assessed.'}

Pitfalls and optional extensions

  • Do not infer stored accuracy from the number of printed digits.
  • Do not use machine epsilon as a universal comparison tolerance.
  • Do not assume every language produces infinity for division by zero; Python raises ZeroDivisionError for 1.0 / 0.0.
  • Try math.nextafter(x, -math.inf) as well as the upward direction, especially at exact powers of two.
  • Repeat the spacing and grouping experiments with another language or numeric type, and document the format before comparing results.
  • Treat the bfloat16 helper as evidence about the declared conversion path, not as evidence about untested mixed-precision arithmetic or named hardware.

Next: Measuring And Comparing Numerical Error.