Learn and teach floating-point rounding#

The problem#

0.1 + 0.2 gives 0.30000000000000004, and students conclude that floating point is approximate and mysterious. To explain it, you need to show what the numbers really are, what the exact answer is, and why the computer picked the result it did.

Why the usual approach falls short#

Printing more digits shows more digits of a rounded number, not the exact value. decimal and fractions can show exact values, but they do not know the formats, roundings and orders that real implementations use.

The solution#

ExactRef shows the exact value of every float, the exact answer of every operation and the rounding that turns one into the other.

What 0.1 really is#

Code
import exactref as er

# The float64 array [0.1]: value((0,)) gives element 0's exact value, and
# to_fraction() writes it as a Fraction
print("exact value:", repr(er.asarray([0.1]).value((0,)).to_fraction()))
Output
exact value: Fraction(3602879701896397, 36028797018963968)

The float nearest 0.1 is a fraction with a power of two below it, slightly more than 0.1.

Why 0.1 + 0.2 gives 0.30000000000000004#

Code
from fractions import Fraction

# ExactRef's reference for 0.1 + 0.2 in float64
s = er.add(er.asarray([0.1]), er.asarray([0.2]))
# The exact sum as a Fraction. Here the enclosure has radius 0, so its
# midpoint is the exact value.
exact = s.exact.value((0,)).midpoint.to_fraction()
print("exact sum:", repr(exact))

# Is the exact sum as far from 0.3 as from 0.30000000000000004? Fraction(x)
# is the exact value of the float x.
print("halfway:", exact - Fraction(0.3) == Fraction(0.30000000000000004) - exact)

# The exact sum rounded once to float64, and Python's own sum
print("rounded:", s.rounded.tolist())
print("Python:", 0.1 + 0.2)
Output
exact sum: Fraction(10808639105689191, 36028797018963968)
halfway: True
rounded: [0.30000000000000004]
Python: 0.30000000000000004

The exact sum of the two floats lies exactly halfway between the floats 0.3 and 0.30000000000000004. Round to nearest, ties to even, picks the one whose last bit is 0, which is 0.30000000000000004. Python is right: that is the correctly rounded sum.

The five rounding modes#

Code
# er.Rounding holds the IEEE 754 rounding modes
for mode in er.Rounding:
    # 1/3 rounded to bfloat16 in this mode
    third = er.asarray([Fraction(1, 3)], "bfloat16", rounding=mode)
    print(f"{mode.name}:", third.tolist())
Output
NEAREST_EVEN: [0.333984375]
NEAREST_AWAY: [0.333984375]
TOWARD_ZERO: [0.33203125]
TOWARD_POSITIVE: [0.333984375]
TOWARD_NEGATIVE: [0.33203125]

How coarse small formats are#

Every float8 e4m3fn value from 1 up to 2, read from its encodings:

Code
# An array of shape (8,) from the encodings 0x38 to 0x3f, read as float8_e4m3fn
values = er.Array.from_encodings("float8_e4m3fn", (8,), list(range(0x38, 0x40)))
print("values:", values.tolist())
Output
values: [1.0, 1.125, 1.25, 1.375, 1.5, 1.625, 1.75, 1.875]

Why the order of additions matters#

Code
# A float64 sum whose exact value is 1
r = er.sum([1e16, 1.0, -1e16])
# The exact sum rounded once to float64
print("rounded:", r.rounded.tolist())
# normal() adds in index order, rounding each partial sum to float64
print("normal:", r.normal().tolist())
Output
rounded: 1.0
normal: 0.0

The exact sum is 1, but adding in index order loses the 1.0, because float64 values near 1e16 are 2 apart. How tolerance works continues from here: why both answers are correct, and how far a correct implementation can be from the exact one.

Subnormals and underflow#

Code
# 1e-45 rounds to the smallest float32 above 0
print("1e-45:", er.asarray([1e-45], "float32").tolist())
# 1e-46 is closer to 0, so it rounds to 0
print("1e-46:", er.asarray([1e-46], "float32").tolist())

# Rounding toward positive infinity instead of to nearest
up = er.Rounding.TOWARD_POSITIVE
print("1e-46 rounded up:", er.asarray([1e-46], "float32", rounding=up).tolist())
Output
1e-45: [1.401298464324817e-45]
1e-46: [0.0]
1e-46 rounded up: [1.401298464324817e-45]