"""Proof: rounding twice can differ from rounding once. Take the
exact value 1 + 2^-24 + 2^-60. Rounded straight to float32 it goes
UP to 1 + 2^-23, because it sits just above the halfway point.
Rounded first to float64, the tiny 2^-60 tail is lost and the
float64 result lands exactly ON the halfway point; rounding that
to float32 then goes DOWN to 1.0 by ties-to-even. Two legal
paths, 2 different answers: double rounding, the failure mode of
every system that computes wider first and stores narrower.
"""
from fractions import Fraction as F
import numpy as np

exact = F(1) + F(1, 2 ** 24) + F(1, 2 ** 60)

# path 1: round once, straight to float32 (decided exactly)
lo, hi = F(1), F(1) + F(1, 2 ** 23)          # float32 neighbors
assert lo < exact < hi
once = np.float32(1.0) if (exact - lo) <= (hi - exact) else \
       np.float32(1.0 + 2.0 ** -23)

# path 2: round to float64 first, then to float32
via64 = float(exact)                          # correctly rounded
twice = np.float32(via64)

print(f"exact value          : 1 + 2^-24 + 2^-60")
print(f"rounded once to fp32 : {once!r}")
print(f"fp64 first           : {via64!r}")
print(f"then to fp32         : {twice!r}")
assert once == np.float32(1.0 + 2.0 ** -23)
assert twice == np.float32(1.0)
print("\nthe 2 paths disagree. the old x87 unit computed in 80-bit")
print("registers and rounded again on every spill to memory, so a")
print("program's bits depended on register allocation; SSE and its")
print("one-width successors retired the problem, and this proof is")
print("why narrow-then-narrower casts still deserve suspicion.")