IEEE 754 Floating Point, in depth: the encoding, correct rounding, special values, ulp comparison, compensated summation and float32 to bfloat16 conversion

The format from first principles

A binary floating-point number is a sign, an exponent and a significand: the value is (-1)^s x 1.f x 2^(e - bias) for normal numbers. Because a normal significand always starts with a 1 bit, that bit is not stored; the fraction field f holds only what follows the binary point. The exponent is stored with a bias so that the field is an unsigned integer, which has a useful consequence: for non-negative values, ordering the bit patterns as integers orders the numbers.

bfloat16 is a truncated binary32: same exponent range, far fewer fraction bits. It follows IEEE conventions but is not one of the standard's formats, which is why conversions to it are defined by libraries and hardware rather than by the standard.

Diagram: one value, bit by bit

Decoding bits exactly

The fastest way to stop guessing is to look at the bits. The function below decodes a binary64 value into an exact rational using Python's Fraction, so you can see precisely which number a literal became.

import functools
import math
import operator
import struct
from fractions import Fraction

def bits64(x: float) -> int:
    return struct.unpack(">Q", struct.pack(">d", x))[0]

def decode64(x: float):
    u = bits64(x)
    sign = u >> 63
    exp = (u >> 52) & 0x7FF
    frac = u & ((1 << 52) - 1)
    if exp == 0x7FF:
        return "nan" if frac else ("-inf" if sign else "inf")
    if exp == 0:                                   # zero or subnormal: no hidden bit
        value = Fraction(frac, 1 << 52) * Fraction(2) ** -1022
    else:                                          # normal: hidden leading 1
        value = (1 + Fraction(frac, 1 << 52)) * Fraction(2) ** (exp - 1023)
    return -value if sign else value

print(hex(bits64(0.1)))            # 0x3fb999999999999a
print(float(decode64(0.1)) == 0.1, decode64(0.1) == Fraction(1, 10))   # True False
print(math.ulp(1.0) == 2.0 ** -52) # True

The literal 0.1 cannot be represented, because one tenth has an infinite binary expansion. The parser picks the nearest representable value, which is slightly above one tenth. Every later surprise, such as 0.1 + 0.2 == 0.30000000000000004, starts from that first rounding.

Rounding: the one guarantee

The core guarantee of IEEE 754 is correct rounding: for addition, subtraction, multiplication, division, square root and fused multiply-add, the result is the exact mathematical result rounded once to the destination format. The default rounding mode is round to nearest, ties to even: if the exact result lies halfway between two neighbours, pick the one whose last bit is zero. Ties to even avoids the steady upward drift that round-half-up would add to long computations.

The distance between adjacent values near x is one ulp (unit in the last place). Near 1.0 in binary64 it is 2^-52, about 2.2e-16, and the worst relative error of one rounding is half of that, 2^-53, the unit roundoff. The guarantee is per operation. A chain of operations accumulates errors, and the order of operations changes the result: floating-point addition is commutative but not associative.

>>> (1e16 + 1.0) - 1e16       # the 1.0 is below half an ulp of 1e16, so it vanishes
0.0
>>> (0.1 + 0.2) + 0.3, 0.1 + (0.2 + 0.3)
(0.6000000000000001, 0.6)      # different order, different rounding, different answer
>>> functools.reduce(operator.add, [0.1] * 10)   # plain left-to-right adds
0.9999999999999999
>>> import functools, operator
>>> import numpy as np
>>> np.float32(2**24 + 1) == np.float32(2**24)   # binary32 has 24 bits of precision
True

Fused multiply-add computes a*b + c with a single rounding instead of two. It makes dot products and polynomial evaluation more accurate, and it also means the same source code can give different bits depending on whether the compiler fused an operation. That is a reproducibility issue, not a bug, and it is why numerical libraries document their FMA use.

Zeros, subnormals, infinities and NaN

• Signed zero. +0.0 and -0.0 compare equal but differ in bits; 1/-0.0 is negative infinity in languages that follow IEEE semantics. Hashing raw bits treats them as different keys.
• Subnormals. When the exponent field is zero, the hidden bit is 0 and the exponent is fixed at its minimum, so values fade gradually toward zero instead of jumping. The smallest positive binary64 subnormal is about 4.9e-324; the smallest normal is about 2.2e-308. Some hardware handles subnormals slowly, and many GPU and SIMD paths offer flush-to-zero modes that trade this gradual underflow for speed.
• Infinities. Overflow rounds to infinity under the default mode, and infinity follows arithmetic rules: inf - inf and 0 * inf are invalid and produce NaN.
• NaN. Not a Number, with a payload in the fraction bits. A NaN compares unequal to everything, including itself, so x != x is the classic test. Quiet NaNs propagate through arithmetic; signaling NaNs raise the invalid exception when used. One NaN in a reduction poisons the whole result, which is useful for detection and painful for debugging.
• Ordering. Ordinary comparisons are a partial order because of NaN, so sorting arrays that contain NaN with comparison-based sorts can give inconsistent results. The standard defines a totalOrder predicate for exactly this case.

Comparing floats

Never test computed floats with == unless the values are exactly representable by construction. Use a tolerance with both a relative part, for large values, and an absolute part, for values near zero where relative error is meaningless. Python's math.isclose does exactly this. For tests of numerical kernels, a distance in ulps is often more informative, because it says how many representable values separate the two answers regardless of magnitude.

def ordered(x: float) -> int:
    """Map a double to an integer whose order matches the float order (and -0.0 == +0.0)."""
    i = struct.unpack("<q", struct.pack("<d", x))[0]
    return i if i >= 0 else -(i & 0x7FFFFFFFFFFFFFFF)

def ulp_distance(a: float, b: float) -> int:
    if math.isnan(a) or math.isnan(b):
        raise ValueError("ulp distance is undefined for NaN")
    return abs(ordered(a) - ordered(b))

assert ulp_distance(1.0, math.nextafter(1.0, 2.0)) == 1
assert ulp_distance(-0.0, 0.0) == 0
assert ulp_distance(0.1 + 0.2, 0.3) == 1

The trick works because of the biased exponent: positive values already sort as integers, and flipping negative values around zero makes the whole line monotonic. A kernel test that asserts ulp_distance(got, want) <= 4 documents its accuracy far better than abs(got - want) < 1e-6.

Summation without losing bits

Summation is where rounding errors pile up. Adding n numbers left to right can have an error bound that grows linearly with n, and adding a small value to a large running sum loses the small value's low bits. Three fixes, in increasing cost: sum in a wider type, sum pairwise (NumPy does this for float arrays, giving error growth closer to log n), or carry the lost bits in a compensation term.

def neumaier_sum(xs):
    """Compensated summation (Kahan-Babuska / Neumaier)."""
    s = 0.0
    comp = 0.0
    for x in xs:
        t = s + x
        if abs(s) >= abs(x):
            comp += (s - t) + x      # low bits of x lost in s + x
        else:
            comp += (x - t) + s      # low bits of s lost in s + x
        s = t
    return s + comp

data = [1.0, 1e100, 1.0, -1e100]
naive = functools.reduce(operator.add, data)
print(naive, neumaier_sum(data), math.fsum(data))   # 0.0 2.0 2.0

Note that since Python 3.12 the built-in sum of floats already uses this compensated algorithm, which is why the example spells out the naive loop; NumPy, C and most GPU kernels do not.

Classic Kahan summation handles a running sum that dominates each term but fails on the example above, where a term dominates the sum; Neumaier's variant handles both. math.fsum returns the correctly rounded exact sum and is the reference to test against. The same reasoning explains why mixed-precision training accumulates matrix products in float32 even when inputs are 16-bit; see mixed precision training.

Variance is the other common trap. The textbook formula E[x^2] - E[x]^2 subtracts two large, nearly equal numbers and can return a negative variance. Welford's one-pass update keeps a running mean and a running sum of squared deviations and stays accurate.

Worked example: float32 to bfloat16

Worked example: converting float32 to bfloat16. Since bfloat16 is the top 16 bits of a float32, the naive conversion truncates, which always rounds toward zero and biases every value. Correct conversion rounds to nearest, ties to even, by adding a bias before the shift. NaN must be handled first: adding the bias to a NaN whose payload sits only in the low bits would carry into the exponent and turn it into infinity.

def f32_to_bf16_bits(x: float) -> int:
    u = struct.unpack("<I", struct.pack("<f", x))[0]
    if (u & 0x7F800000) == 0x7F800000 and (u & 0x007FFFFF):
        return ((u >> 16) | 0x0040) & 0xFFFF       # keep sign, force a quiet NaN
    lsb = (u >> 16) & 1                            # last bit that survives
    u += 0x7FFF + lsb                              # ties go to the even neighbour
    return (u >> 16) & 0xFFFF

def bf16_bits_to_f32(b: int) -> float:
    return struct.unpack("<f", struct.pack("<I", b << 16))[0]

for v in (1.0, 0.1, 3.0e38, -2.5):
    b = f32_to_bf16_bits(v)
    print(v, hex(b), bf16_bits_to_f32(b))
# 0.1 -> 0x3dcd -> 0.10009765625 ; 3.0e38 -> 0x7f62 (still finite)

Trace 0.1: its float32 bits are 0x3DCCCCCD. The low 16 bits 0xCCCD exceed half of 0x10000, so the result rounds up from 0x3DCC to 0x3DCD, which decodes to 0.10009765625. Truncation would have given 0x3DCC, 0.099609375, which is further away. Values above the largest bfloat16 that round up overflow correctly to infinity, because the carry lands in the exponent field exactly as the format intends.

Failure modes

• Equality on computed values. Loops that stop at t == 1.0 after adding 0.1 ten times never stop. Count steps with integers and derive the float.
• Catastrophic cancellation. Subtracting nearly equal values leaves only the rounding noise. Rearrange formulas (Welford, log1p, expm1, hypot) instead of adding precision.
• Money in binary floats. Decimal fractions are not representable. Use integer minor units or a decimal type.
• Non-reproducible reductions. Parallel sums reorder additions, so results change with thread count. Fix the reduction order or use compensated or exact sums where bitwise reproducibility matters.
• Fast-math flags. Compiler options that assume no NaN or infinity can delete your x != x checks.
• Silent NaN spread. One NaN in a loss or a metric aggregate erases the result. Check for NaN at boundaries and log the first occurrence with its inputs.

Trade-offs

What to do next

1. Decode a few of your own constants with the decode function above and note which are inexact.
2. Replace float equality in tests with math.isclose or a ulp-distance assertion that states the expected accuracy.
3. Find your largest reductions and check their accumulator type and order; test them against math.fsum.
4. Audit variance and difference computations for cancellation and switch to stable formulas.
5. If you convert to 16-bit formats, confirm the converter rounds to nearest even and handles NaN, then read bfloat16 in depth and FP8 formats.