IEEE 754 anatomy, ulp, round to nearest even, and bf16 and fp16 emulation

Overview

Module	M09.1 · build · Python · Pass 2 · 2 to 3 h
You build	python/tinyllm/num/fp.py: decompose_f32, compose_f32, ulp, f32_to_bf16_bits, bf16_bits_to_f32, round_to_bf16, round_to_fp16
Contract	course/contracts/py/tinyllm/num/fp.pyi
Tests	course/tests/M09.1/test_fp.py (what they check: section 4)
Needs	nothing to build first. Reading: M00.1 exponents and logarithms (powers of two, $\log_2$)
Used by	L0.6 reads and writes BF16 safetensors through f32_to_bf16_bits and bf16_bits_to_f32 · later L7.9 loads bf16 checkpoints through bf16_bits_to_f32, L11.1 rounds weights and activations with round_to_bf16, M09.4 builds fp8 and MX formats on the same bit arithmetic
Milestone	MS-P2 (the foundations gate)
Optional depth	Goldberg, “What Every Computer Scientist Should Know About Floating-Point Arithmetic” (ACM Computing Surveys, 1991), sections 1 and 2; Muller et al., Handbook of Floating-Point Arithmetic (2nd ed.), ch. 2 and 3; IEEE Std 754-2019, section 3

Key Takeaways

• A float32 is three integers: sign $s$, biased exponent $E$, and mantissa $M$; a normal value is $(-1)^s (1 + M/2^{23})\, 2^{E - 127}$, and $E = 0$ and $E = 255$ hold the zeros, subnormals, infinities, and NaNs (test_decompose_hand_example, test_decompose_every_class).
• The gap between neighbouring floats, one ulp, is $2^{e - p + 1}$: it doubles at every power of two and stops shrinking in the subnormal range (test_ulp_matches_numpy_spacing).
• bfloat16 is the top half of a float32, so rounding to it is integer arithmetic on the bits: add $\mathtt{0x7FFF}$ plus the lowest kept bit, then shift. Ties go to the even neighbour, and a carry walks into the exponent exactly as it should (test_bf16_matches_exact_rational_golden).
• The trick is wrong for NaN, which can turn into infinity; every NaN must map to the quiet code 0x7FC0 (test_bf16_nan_stays_nan).
• Codes are decoded with a bit view, never a value cast: astype turns the code 0x3F80 into $16256.0$ instead of $1.0$ (test_bf16_all_codes_roundtrip).

How to work this chapter

ol start M09.1              # stubs python/tinyllm/num/fp.py into your repo
ol tests M09.1              # read the test catalog first: rung R0, you write no tests here
ol check M09.1              # exit code is the verdict
ol diff  M09.1              # after passing: your code against the reference

---

1. Why now

Everything your tracer stores is float32: L0.0 writes the bigram’s weights as F32 in safetensors, and M03.1 multiplies them in C. Pass 2 starts to lean on what float32 actually is. The stable softmax of M09.2 exists because float32 overflows at $e^{88.72}$, and that number comes from the exponent field. Later passes go below 32 bits. The SmolLM2 checkpoint L7.9 loads is stored as BF16, two bytes per weight, and the obvious loader, np.frombuffer(raw, np.uint16).astype(np.float32), gives you a model whose weight 1.0 reads as 16256.0. Its first forward pass returns NaN and nothing crashes. Mixed-precision training (L11.1) rounds every weight to bf16, and the obvious rounding (keep the top 16 bits) is biased: it always rounds toward zero, and over thousands of steps the weights drift. This module takes the float32 format apart bit by bit and rebuilds the two 16-bit formats from it, so the loader and the rounding are right the first time.

2. Principles

Symbol	Meaning	Type
$s$	sign bit: 0 for positive, 1 for negative	int, 0 or 1
$E$	biased exponent field: 8 bits in float32	int, 0 to 255
$M$	mantissa (fraction) field: 23 bits in float32	int, 0 to $2^{23} - 1$
$e = E - 127$	the unbiased exponent of a normal float32	int, -126 to 127
$p$	precision: significand bits including the implicit leading 1	24 (f32), 11 (f16), 8 (bf16)
$e_{\min}$	smallest normal exponent	-126 (f32, bf16), -14 (f16)
$\mathrm{ulp}(x)$	the gap from $\lvert x \rvert$ to the next larger magnitude in the format	float64
$u = 2^{-p}$	unit roundoff: the largest relative error of one rounding	float
$\mathrm{rn}(x)$	$x$ rounded to the nearest representable value, ties to even	a value of the format
$b$	the 32 bits of a float32 read as an unsigned integer (a uint32 view)	int

2.1 The binary32 encoding

A float32 occupies 32 bits: bit 31 is the sign $s$, bits 30 to 23 the biased exponent $E$, bits 22 to 0 the mantissa $M$. Reading the same 32 bits as an unsigned integer gives $b = s \cdot 2^{31} + E \cdot 2^{23} + M$, so the fields are b >> 31, (b >> 23) & 0xFF, and b & 0x7FFFFF. The value depends on $E$:

$E$	class	value
1 to 254	normal	$(-1)^s \left(1 + \frac{M}{2^{23}}\right) 2^{E - 127}$
0	zero ($M = 0$) or subnormal ($M \ne 0$)	$(-1)^s \frac{M}{2^{23}}\, 2^{-126}$
255	infinity ($M = 0$) or NaN ($M \ne 0$)	$\pm\infty$, or not a number

A normal number has $p = 24$ significant bits: the 23 stored ones plus a leading 1 that is implied, not stored. The bias 127 lets the exponent field be an unsigned integer, which has a useful consequence: for non-negative floats, ordering the bit patterns as integers orders the values. The largest finite float32 is $E = 254$, $M = 2^{23} - 1$, which is $(2 - 2^{-23})\, 2^{127} \approx 3.40 \times 10^{38} = e^{88.72}$. The smallest normal is $2^{-126}$. Below it, $E = 0$ drops the implicit 1 and keeps the exponent at $-126$, so the subnormals fill the gap down to $2^{-149}$ in equal steps. There are two zeros, $+0$ and $-0$, which compare equal. A NaN is any pattern with $E = 255$ and $M \ne 0$; there are about $2^{24}$ of them.

2.2 Spacing and the ulp

Between $2^e$ and $2^{e+1}$ (one binade) a format with precision $p$ has $2^{p-1}$ equally spaced values, so their spacing is

$\mathrm{ulp}(x) = 2^{\,e - p + 1}, \qquad e = \max(\lfloor \log_2 \lvert x \rvert \rfloor,\ e_{\min}).$

The spacing doubles at every power of two and is constant across the subnormal range, which is why $e$ is clamped at $e_{\min}$. At $x = 1$: $2^{-23}$ in float32, $2^{-10}$ in float16, $2^{-7}$ in bfloat16. At $x = 0$ the gap to the next value is the smallest subnormal, $2^{e_{\min} - p + 1}$. Infinity and NaN have no next value, so their ulp is NaN.

Rounding to nearest moves a value by at most half an ulp, and within a binade the ulp is at most $2^{-(p-1)} \lvert x \rvert$, so one rounding has a relative error of at most $u = 2^{-p}$: $\mathrm{rn}(x) = x(1 + \delta)$ with $\lvert \delta \rvert \le u$, for any $x$ in the normal range. That one fact is the model of every error bound in this track (M09.2, M09.3).

Computing $\lfloor \log_2 \lvert x \rvert \rfloor$ needs care. np.log2(x) returns a rounded float, and for $x$ a hair below $2^k$ the true value $k - \varepsilon$ rounds up to exactly $k$, giving the next binade’s ulp. The exponent is already stored exactly: np.frexp(x) returns $(m, k)$ with $x = m\, 2^k$ and $m \in [0.5, 1)$, so $\lfloor \log_2 \lvert x \rvert \rfloor = k - 1$ with no rounding at all.

2.3 Round to nearest even, on the bits, for bfloat16

bfloat16 has the float32 layout cut in half: 1 sign bit, the same 8 exponent bits with the same bias, and 7 mantissa bits ($p = 8$). Every bf16 code is the top 16 bits of a float32, and decoding is a shift: bf16_bits_to_f32(c) is the float32 whose bits are c << 16. Because the range is float32’s, nothing a float32 model produces overflows on the way down; only precision is lost.

Encoding has to round. IEEE 754’s default is round to nearest, ties to even: take the nearest representable value, and when the value is exactly halfway between two, take the one whose last kept bit is 0. Ties to even is unbiased: half of all ties round up and half round down, where “ties away from zero” pushes every tie the same way.

Split the float32 bits into the half we keep and the half we drop: $h = b \gg 16$ and $\ell = b \mathbin{\&} \mathtt{0xFFFF}$. The halfway point of the dropped half is $\ell = \mathtt{0x8000}$. The rule:

dropped half $\ell$	result
below 0x8000	$h$ (round down)
above 0x8000	$h + 1$ (round up)
exactly 0x8000	$h$ if $h$ is even, $h + 1$ if odd

All three rows are one integer expression:

$\mathrm{bf16}(b) = (b + \mathtt{0x7FFF} + (h \mathbin{\&} 1)) \gg 16.$

When $\ell < \mathtt{0x8000}$, adding at most $\mathtt{0x8000}$ keeps the low half below $\mathtt{0x10000}$, so no carry reaches $h$. When $\ell > \mathtt{0x8000}$, adding $\mathtt{0x7FFF}$ already carries. When $\ell = \mathtt{0x8000}$, the sum is $\mathtt{0xFFFF}$ plus the parity bit, so it carries exactly when $h$ is odd.

The carry needs no special case. If the 7 kept mantissa bits are all ones, adding 1 to $h$ clears them and increments the exponent field: the next code is the first value of the next binade, which is the correct rounded result. At the top, the largest finite float32 0x7F7FFFFF rounds to 0x7F80, which is $+\infty$: exactly IEEE’s rule, since $(2 - 2^{-23})\, 2^{127}$ is nearer to $2^{128}$ than to the largest bf16 value. Infinity itself has $\ell = 0$ and stays infinity. Work in uint64 (or check for it) so that the addition cannot wrap at 0xFFFFFFFF.

NaN breaks it. A NaN may have its only set mantissa bits in the dropped half: 0x7F800001 has $h = \mathtt{0x7F80}$ and $\ell = 1$, and so does every NaN whose payload is small. Truncated, it becomes 0x7F80, infinity. Rounded, 0x7F800001 + 0x7FFF = 0x7F808000, still infinity. A NaN loss that turns into infinity no longer trips np.isnan. So NaNs are mapped first, to the single quiet NaN code 0x7FC0, matching the C conversions in tinyllm/numerics.h.

2.4 binary16: a different exponent range

float16 (IEEE binary16) has 1 sign bit, 5 exponent bits with bias 15, and 10 mantissa bits ($p = 11$, $e_{\min} = -14$). Its range is tiny: the largest finite value is $(2 - 2^{-10})\, 2^{15} = 65504$, and the smallest subnormal is $2^{-24}$. So a float32 can overflow it, land in its subnormal range, or underflow it entirely, and the exponent field has to be re-biased. The top-half trick does not apply.

Write a finite $\lvert x \rvert = \sigma \cdot 2^{E' - 23}$ with an integer significand $\sigma < 2^{24}$ ($\sigma = 2^{23} + M$ and $E' = E - 127$ for a normal float32, $\sigma = M$ and $E' = -126$ for a subnormal). binary16 spaces its values $2^{E_h - 10}$ apart with $E_h = \max(E', -14)$. Dividing $\lvert x \rvert$ by that spacing gives $q = \sigma / 2^{t}$ with $t = 13 + (E_h - E') \ge 13$: shift right by $t$ and round the remainder to nearest even, exactly as for bf16 but with a variable cut. Then

$\text{code} = ((E_h + 14) \ll 10) + q.$

For a normal result ($2^{10} \le q < 2^{11}$) this equals the field form $((E_h + 15) \ll 10) \mid (q - 2^{10})$; for a subnormal ($E_h = -14$, $q < 2^{10}$) it is just $q$; and a $q$ that rounded up to $2^{11}$ carries into the exponent on its own, as in bf16. Any code at or above 0x7C00 is infinity (overflow), and NaN gets 0x7E00. numpy’s float32.astype(float16) implements the same IEEE rounding, which gives the tests an independent oracle.

2.5 Bit views versus value casts

np.asarray(x, np.float32).view(np.uint32) reinterprets the same bytes as integers and changes no bit; .astype(np.uint32) converts the value (1.0 becomes 1). Every function here moves between the two worlds with views. Decoding bf16 is the case people get wrong: the codes arrive as uint16, and the value of code $c$ is the float32 with bits $c \ll 16$, not the number $c$.

3. Worked example by hand

Anatomy of $-6.25$. $6.25 = 110.01_2 = 1.1001_2 \times 2^2$. So $s = 1$, $E = 2 + 127 = 129 = 1000\,0001_2$, and the mantissa is the fraction after the leading 1, $.1001_2$, padded to 23 bits: $100\,1000\,0000\,0000\,0000\,0000_2 = \mathtt{0x480000}$. All 32 bits: 1 10000001 10010000000000000000000 = 0xC0C80000. decompose_f32(-6.25) == (1, 129, 0x480000). Check: $(1 + \mathtt{0x480000}/2^{23}) \cdot 2^{2} = (1 + 0.5625) \cdot 4 = 6.25$.

ulp. $6.25$ lies in the binade $[4, 8)$, $e = 2$, so $\mathrm{ulp}(6.25) = 2^{2 - 23} = 2^{-21}$ in float32, $2^{2 - 10} = 2^{-8}$ in float16, and $2^{2 - 7} = 2^{-5}$ in bfloat16. At 1.0 the three are $2^{-23}$, $2^{-10}$, $2^{-7}$ (test_ulp_hand_values).

$0.1$ to bfloat16. $0.1$ has no finite binary expansion. float32(0.1) has bits 0x3DCCCCCD and the value $0.100000001490116119384765625$. The kept half is $h = \mathtt{0x3DCC}$, the dropped half $\ell = \mathtt{0xCCCD}$:

step	value
$\ell$ vs 0x8000	0xCCCD > 0x8000: round up
$h \mathbin{\&} 1$	0 (0x3DCC is even)
$b + \mathtt{0x7FFF} + 0$	0x3DCCCCCD + 0x7FFF = 0x3DCD4CCC
$\gg 16$	0x3DCD

Decode 0x3DCD = 0 01111011 1001101: $E = 123$, $e = -4$, mantissa $1001101_2 = 77$, value $(1 + 77/128)\, 2^{-4} = 205/2048 = 0.10009765625$. The error, about $9.8 \times 10^{-5}$, is below half the bf16 ulp at $0.1$, $2^{-12} \approx 2.4 \times 10^{-4}$.

Ties. $1 + 2^{-8}$ has bits 0x3F808000: $h = \mathtt{0x3F80}$ (that is 1.0), $\ell = \mathtt{0x8000}$, exactly halfway to $1 + 2^{-7}$. $h$ is even, so the sum 0x3F808000 + 0x7FFF = 0x3F80FFFF does not carry: the result is 0x3F80, $1.0$. Next, $1 + 3 \cdot 2^{-8}$ has bits 0x3F818000: halfway between 0x3F81 (odd) and 0x3F82 (even). Now the parity bit is 1, 0x3F818000 + 0x7FFF + 1 = 0x3F820000, and the result is 0x3F82 $= 1 + 2^{-6}$. One tie went down and one went up, both to the even code (test_bf16_hand_examples).

$0.1$ to float16. $e = -4 \ge -14$, so the result is normal and the cut is at $t = 13$. The float32 mantissa 0x4CCCCD is $100\,1100\,1100\,1100\,1100\,1101_2$; the top 10 bits are $10\,0110\,0110_2 = \mathtt{0x266}$ and the 13 dropped bits are $0\,1100\,1100\,1101_2 = \mathtt{0x0CCD}$, below the halfway point 0x1000, so it rounds down. The code is $((-4 + 15) \ll 10) \mid \mathtt{0x266} = \mathtt{0x2C00} \mid \mathtt{0x266} = \mathtt{0x2E66}$, the value $(1 + 614/1024)\, 2^{-4} = 0.0999755859375$. float16 lands closer to $0.1$ than bf16 did: 3 more mantissa bits.

Edges of float16. $65519$ is below the midpoint $65520$ between $65504$ and the would-be next value $65536$, so it rounds to $65504$; $65520$ is the tie, and ties to even picks $65536$, which does not exist, so the result is $+\infty$. At the bottom, $2^{-25}$ is half the smallest subnormal $2^{-24}$: a tie between $0$ and $2^{-24}$, and the even choice is $0$ (test_fp16_edges).

4. The interface

def decompose_f32(x: float) -> tuple[int, int, int]: ...          # (sign, biased exponent, mantissa)
def compose_f32(sign: int, exponent: int, mantissa: int) -> float: ...
def ulp(x: ArrayLike, dtype: Literal["f32", "f16", "bf16"]) -> NDArray: ...   # float64
def f32_to_bf16_bits(x: ArrayLike) -> NDArray: ...                # uint16, RNE, NaN -> 0x7FC0
def bf16_bits_to_f32(u16: ArrayLike) -> NDArray: ...              # float32, exact
def round_to_bf16(x: ArrayLike) -> NDArray: ...                   # float32 in, float32 out
def round_to_fp16(x: ArrayLike) -> NDArray: ...                   # float32, same bits as numpy's cast

Array functions accept anything numpy converts, convert it to float32 first, and keep the input’s shape (a scalar gives a 0-d array). L7.9 reads a BF16 tensor as bf16_bits_to_f32(np.frombuffer(raw, "<u2").reshape(shape)); L11.1 keeps float32 master weights and computes with round_to_bf16(w).

What the tests check

Test	KIND	Checks	Why it matters downstream
test_decompose_hand_example	unit, smoke	section 3’s $-6.25 \to (1, 129, \mathtt{0x480000})$, plus 1.0 and $-2.5$	you and the test agree on the fields
test_decompose_every_class	boundary	$\pm 0$, the smallest subnormal and normal, the largest finite, $\infty$, NaN	a checkpoint can hold any of them
test_compose_roundtrip	property	compose inverts decompose on 3000 floats; out-of-range fields raise	the three fields describe the number completely
test_ulp_hand_values	unit, smoke	$\mathrm{ulp}(1)$ in each format, $\mathrm{ulp}(0)$, NaN for $\pm\infty$ and NaN, unknown dtype	the relative precision of each format
test_ulp_matches_numpy_spacing	differential	float32 and float16 ulps equal np.spacing, subnormals included	the spacing M09.3 builds tolerances from
test_ulp_just_below_a_power_of_two	boundary	nextafter(1024, 0) and nextafter(2**100, 0) stay in the lower binade	np.log2 rounds up there
test_bf16_hand_examples	unit, smoke	section 3: 0.1 -> 0x3DCD, the two ties	the rounding rule itself
test_bf16_matches_exact_rational_golden	golden	10 156 float32 patterns against rounding done with exact fractions	every binade, tie, and overflow
test_bf16_nan_stays_nan	boundary	four NaNs, including 0x7F800001, give 0x7FC0	a NaN loss stays visible in L11.1
test_bf16_all_codes_roundtrip	property	all 65536 codes decode exactly and encode back	L7.9 loads every weight exactly
test_rounding_is_idempotent_and_monotone	property	rounding twice changes nothing; sorted inputs stay sorted	clipping and comparisons after rounding
test_fp16_matches_numpy_bitwise	differential, golden	the golden inputs and $10^5$ random patterns equal numpy’s cast bit for bit	an independent IEEE implementation agrees
test_fp16_edges	boundary	$65504$, $65519 \to 65504$, $65520 \to \infty$, subnormals, $2^{-25} \to 0$	why fp16 training needs loss scaling and bf16 does not
test_shapes_and_scalars	boundary	matrices keep their shape, a scalar gives a 0-d array, the sign survives	weight matrices go through unchanged

The golden file course/fixtures/M09.1/round_f32.npz was written by course/oracle/M09.1/lowp_golden.py, which rounds each float32 value as an exact fractions.Fraction and checks its float16 column against numpy.

5. Pitfalls

Pitfall	Symptom	Caught by
keeping the top 16 bits (truncation)	every value rounds toward zero; 0.1 gives 0x3DCC	test_bf16_hand_examples (mutant s01)
adding 0x8000 (ties away from zero)	1 + 2**-8 rounds to 0x3F81, a biased tie	test_bf16_hand_examples (mutant s02)
no special case for NaN	0x7F800001 becomes 0x7F80, infinity	test_bf16_nan_stays_nan (mutant s03)
$\lfloor \log_2 x \rfloor$ through np.log2	the ulp just below a power of two is twice too large	test_ulp_just_below_a_power_of_two (mutant s04)
no clamp at $e_{\min}$	subnormal ulps keep shrinking past the real spacing	test_ulp_matches_numpy_spacing (mutant s05)
flushing float16 subnormals to zero	values below $2^{-14}$ vanish	test_fp16_edges (mutant s06)
no clamp at 0x7C00	values past 65504 wrap into NaN codes	test_fp16_edges (mutant s07)
float16 ties rounded up	half of all ties off by one ulp	test_fp16_matches_numpy_bitwise (mutant s08)
astype instead of a bit view	code 0x3F80 reads as $16256.0$	test_bf16_all_codes_roundtrip (mutant s09)
exponent not masked with 0xFF	the sign bit leaks into the exponent of negative numbers	test_decompose_hand_example (mutant s10)

6. Where it’s used next

Direction	Module	How it uses this
Back	M00.1	powers of two and $\log_2$ (reading)
Forward	M09.2	where float32 overflows and underflows, and the unit roundoff behind its summation bounds (reading)
Forward	L0.6	safetensors for every dtype: BF16 tensors are encoded and decoded with this module’s bit converters
Forward	M09.4	fp8 E4M3 and E5M2 and MX block formats: the same bit-level rounding with fewer bits, in Python and C
Forward	L7.9	loads the SmolLM2 BF16 safetensors through bf16_bits_to_f32
Forward	L11.1	mixed precision: float32 master weights, round_to_bf16 for the compute copy

If you skip this module, ol check L0.6 stops with needs M09.1: build it, or pass --ref-deps.

Going further

Your piece	Production equivalent	What it adds	Where to look
f32_to_bf16_bits	PyTorch c10::BFloat16	the same add-0x7FFF-plus-parity rounding, in C++ and in CUDA for every kernel	c10/util/BFloat16.h, round_to_nearest_even
round_to_fp16	the F16C instructions (vcvtps2ph), CUDA __float2half_rn	the conversion in one hardware instruction, with the rounding mode as an operand	Intel SDM, VCVTPS2PH; CUDA Math API, half precision intrinsics
round_to_bf16 in training	mixed-precision training with loss scaling	float32 master weights, low-precision compute, dynamic loss scaling for fp16’s narrow range	Micikevicius et al., “Mixed Precision Training” (ICLR 2018)
bf16 and fp16	OCP FP8 and MX formats	8-bit and block-scaled formats built on the same encoding rules	OCP Microscaling Formats (MX) v1.0 spec; M09.4