Machine numbers¶

Level: 201 → 301 · for anyone who has used a float

One line: A float is not a blurry real number but an exact member of a finite set picked by four numbers — radix, precision and two exponent limits — and every operation is exact arithmetic followed by a rounding function, which decides, law by law, what survives of real arithmetic.

A finite set standing in for the real line¶

The first three lessons of this chapter were about numbers that came from a measurement, whose later digits were never known. This one is about a different gap: the one between mathematics and the machine.

A computer can hold only finitely many different numbers, and the real line has uncountably many points. So almost all of them are missing, and not by bad luck. A finite set has measure zero, so a real number picked at random is a machine number with probability zero.

What is striking is how little it takes to describe the set the machine keeps instead. Four numbers do it:

symbol	name	binary64, Python's `float`
β	radix	2
p	precision, in digits	53
e_min	smallest exponent	−1022
e_max	largest exponent	1023

A normal member is ±d₀.d₁d₂…d₍ₚ₋₁₎ × β^e, with d₀ ≠ 0 and e between e_min and e_max: scientific notation with exactly p significant digits and a limited exponent. That is the whole definition.

Every member is an exact number. The toy member 0.625 below is precisely five eighths. The approximation is not inside the members; it is in the step that forces a real result back into the set.

The Handbook of Floating-Point Arithmetic opens by making the same point: a finite set cannot keep every property of real arithmetic, so every format trades speed, accuracy, range, simplicity and memory against each other. This page takes that literally. Section 4 finds out which properties are kept, and section 7 what each of the four numbers buys.

The program builds a system small enough to print whole — β = 2, p = 3, e from −1 to 2 — tests the laws of arithmetic against every one of its members, and then shows that Python's float is the same construction with bigger numbers.

What the program prints¶

Verified output of machine_numbers.py — regenerated by tools/run_examples.py, never hand-typed.

1. FOUR NUMBERS PICK A SET
   A toy system, small enough to print whole:
     radix beta = 2    precision p = 3 digits    exponent e from -1 to 2
   A member is  d0.d1d2 (in binary) x 2^e,  with d0 = 1:

     significand     e = -1    e = 0    e = 1    e = 2
       1.00 = 1         0.5        1        2        4
       1.01 = 1.25    0.625     1.25      2.5        5
       1.10 = 1.5      0.75      1.5        3        6
       1.11 = 1.75    0.875     1.75      3.5        7

   16 positive, 16 negative and zero: 33 members.
   Per sign: (beta - 1) x beta^(p-1) x (number of exponents) = 1 x 4 x 4 = 16.
   Largest 7, smallest positive 0.5.
   Every other real number -- 0.1, 1/3, 0.51, pi -- is simply not in the set.

2. THE GAPS GROW WITH THE NUMBERS
   Every multiple of 1/8 from 0 to 7, marked where a member sits:

     0       1       2       3       4       5       6       7
     |...|||||.|.|.|.|...|...|...|...|.......|.......|.......|

   Between one power of 2 and the next, the members are evenly spaced:
     from 0.5 to 1     gap 0.125
     from   1 to 2     gap 0.25
     from   2 to 4     gap 0.5
     from   4 to 7     gap 1
   The gap doubles when the number doubles, so the gap RELATIVE to the
   number stays between 1/7 and 1/4 everywhere. A float keeps a fixed
   number of significant digits, never a fixed number of decimal places.

   THE HOLE AT ZERO. Below 0.5 there is nothing but 0. So two different
   members can have a difference too small to be a member:
     0.625 (-) 0.5  =  round(0.125)  =  0      although 0.625 != 0.5
   IEEE 754 fills the hole with SUBNORMALS: d0 = 0 is allowed at e = -1.

     |||||||||.|.|.|.|...|...|...|...|.......|.......|.......|

     0.625 (-) 0.5  =  round(0.125)  =  0.125
   From here on the toy system has them: 19 per sign, 39 members.

3. ROUNDING IS A FUNCTION, AND IEEE 754 DEFINES FIVE
   Each one sends every real number to a member, or to an infinity:
     RD toward -inf    RU toward +inf    RZ toward zero
     RN to nearest, a tie going to the EVEN significand or AWAY from zero

     x          RD        RU        RZ        RN_even   RN_away
     1/10       0         0.125     0         0.125     0.125
     1/3        0.25      0.375     0.25      0.375     0.375
     1.625      1.5       1.75      1.5       1.5       1.75
     -1.625     -1.75     -1.5      -1.5      -1.5      -1.75
     7.25       7         +inf      7         7         7
     7.5        7         +inf      7         +inf      +inf
     0.03125    0         0.125     0         0         0

   1.625 is a tie: 1.5 = 1.10 ends in 0 and 1.75 = 1.11 does not.
   7.5 is half a gap past the largest member, so both RNs overflow.

   Checked on every multiple of 1/64 from -8 to 8 (1,025 values):
     RD(x) <= x <= RU(x)                                True
     both RNs always return RD(x) or RU(x)              True
     RZ(x) is RD(x) when x >= 0, and RU(x) when x < 0   True

   Python's float only ever uses RN even. Its decimal module has all five:
     x      ROUND_FLOOR     ROUND_CEILING   ROUND_DOWN      ROUND_HALF_EVEN ROUND_HALF_UP
     2.5    2               3               2               2               3
     -2.5   -3              -2              -2              -2              -3
     2.4    2               3               2               2               2
   (ROUND_HALF_UP is the decimal module's name for ties away from zero.)

   An operation is the exact result, then a rounding function:
     a (+) b = RN(a + b),  and the same for (-), (x) and (/).
     1.25 (+) 0.625 = RN(1.875) = 2      0.125 (x) 0.25 = RN(0.03125) = 0
   From here on, every operation rounds to nearest, ties to even.

4. WHICH LAWS OF ARITHMETIC SURVIVE
   Tested on every member (39), every pair (1,521) and every triple
   (59,319) of the toy system. Cases where a step overflows are skipped.

   a (+) b  =  b (+) a
       holds: 0 failures in 1,401 pairs
   a (x) b  =  b (x) a
       holds: 0 failures in 1,233 pairs
   a <= b   implies   a (+) c <= b (+) c
       holds: 0 failures in 26,162 triples
   a (-) b = 0   implies   a = b
       holds: 0 failures in 1,401 pairs
       ...but without subnormals: 44 failures in 969 pairs, e.g. 0.625 (-) 0.5 = 0
   b/2 <= a <= 2b   implies   a (-) b is exact  (Sterbenz)
       holds: 0 failures in 139 pairs
       ...but without subnormals: 28 failures in 124 pairs

   (a (+) b) (+) c  =  a (+) (b (+) c)
       FAILS: 13,972 failures in 47,823 triples
       (1 (+) 0.125) (+) 0.125 = 1,  but  1 (+) (0.125 (+) 0.125) = 1.25
   a (x) (b (+) c)  =  (a (x) b) (+) (a (x) c)
       FAILS: 9,444 failures in 37,463 triples
       1.5 (x) (1 (+) 0.125) = 1.5,  but  (1.5 (x) 1) (+) (1.5 (x) 0.125) = 1.75
   b != 0   implies   a (+) b != a
       FAILS: 80 failures in 1,362 pairs
       4 (+) 0.125 = 4
   a (x) (1 (/) a)  =  1
       FAILS: 8 failures in 36 members
       7 (x) (1 (/) 7) = 7 (x) 0.125 = 0.875

   The survivors each compare one exact result with itself: a + b and
   b + a are the same real number before any rounding happens. The
   failures compare two ROUTES, and each route rounds in different places.

5. THE FLOAT YOU ALREADY USE IS THE SAME CONSTRUCTION
   Read off this machine:  beta = 2,  p = 53,  e from -1022 to 1023
     members              18,437,736,874,454,810,623
       equal to 2^64 - 2^53 - 1?      True
     largest              1.7976931348623157e+308
       equal to (2 - 2^-52) x 2^1023? True
     smallest normal      2.2250738585072014e-308  = 2^-1022: True
     smallest subnormal   5e-324  = 2^-1074: True
     gap above 1.0        2.220446049250313e-16
     gap above 1e+16      2.0
     gap above 1e+308     1.99584030953472e+292

   The laws of section 4, on real floats:
     0.1 + 0.2 == 0.2 + 0.1       True
     (0.1 + 0.2) + 0.3            0.6000000000000001
     0.1 + (0.2 + 0.3)            0.6
     100.0 * (0.1 + 0.2)          30.000000000000004
     100.0 * 0.1 + 100.0 * 0.2    30.0
     2.0**53 + 1.0 == 2.0**53     True
     49.0 * (1.0 / 49.0)          0.9999999999999999
     (0.1 + 0.2) - 0.3            5.551115123125783e-17
       ...exactly the difference of the two stored values? True

   Which fractions can be members at ANY precision? In radix beta, a
   fraction in lowest terms can be, only if every prime factor of its
   denominator also divides beta:

                radix 2   radix 10  radix 3
     1/10       no        yes       no
     1/3        no        no        yes
     3/8        yes       yes       no
     1/6        no        no        no

6. OFF THE ENDS OF THE RANGE
   The largest float is about 1.8e308. Factorials outrun it quickly:
     170!  as a float  = 7.257415615307994e+306
     171!  as a float  = inf
     1000! has 2,568 digits, 1001! has 2,571

   1001! / 1000! is exactly 1001. With floats:
     1001! -> inf    1000! -> inf    inf / inf = nan
   If overflow stopped at the largest float instead of going to inf:
     largest / largest = 1.0     a quiet, plausible, wrong answer

   Three ways to stay inside the set:
     cancel first      1001! / 1000!  =  1001 x 1000! / 1000!  =  1001
     exact integers    factorial(1001) // factorial(1000)  =  1001
     logarithms        exp(lgamma(1002) - lgamma(1001))  =  1001   to 9 s.f.

   And off the bottom, where the smallest subnormal is 5e-324:
     5e-324 / 2  =  0.0     a tie between 0 and 5e-324, and 0 is even

7. CHOOSING THE FOUR NUMBERS IS THE COMPROMISE
                bits  radix    p    emin   emax   largest     members
     binary16     16      2   11     -14     15   6.55E+4     6.35e+04
     binary32     32      2   24    -126    127   3.40E+38    4.28e+09
     binary64     64      2   53   -1022   1023   1.79E+308   1.84e+19
     decimal64    64     10   16    -383    384   9.99E+384   1.38e+19
   Brute force on binary16: its 65,536 bit patterns give 63,487 distinct
   finite values. The count formula gives 63,487.

   What each one stores for a few numbers, as a relative error:
                1/10        1/3         2^24 + 1    70,000
     binary16   -2.4e-04    -2.4e-04    overflow    overflow
     binary32   +1.5e-08    +3.0e-08    -6.0e-08    exact
     binary64   +5.6e-17    -5.6e-17    exact       exact
     decimal64  exact       -1.0e-16    exact       exact

   Relative gap just above 1, and just below the next power of the radix:
     binary64    2.2e-16  and  1.1e-16    varies 2x
     decimal64   1.0e-15  and  1.0e-16    varies 10x

   More bits buy precision and range, and cost memory and speed. Radix 10
   makes 1/10 exact, but its gaps wobble 10x instead of 2x, and the same
   64 bits give fewer distinct values. Every choice gives something up.

The gaps are a ruler whose marks spread out¶

Section 2's ruler is the picture to keep. Between one power of 2 and the next, the members are evenly spaced, and at each power of 2 the spacing doubles. The absolute gap grows with the numbers while the relative gap stays inside a fixed band: 1/7 to 1/4 in the toy system, about 1.1 × 10⁻¹⁶ to 2.2 × 10⁻¹⁶ in binary64.

That is exactly what significant figures are. A float stores a fixed number of significant digits, like :.17g, and never a fixed number of decimal places, like :.2f. Between 2⁵³ and 2⁵⁴ the gap is 2, so no odd integer there is a member, which is why 2.0**53 + 1.0 == 2.0**53.

The hole at zero is the one place the ruler has to be patched. With only normal members, nothing lies between 0 and the smallest normal number, so two different members can differ by less than that: 0.625 − 0.5 rounds to 0 although the two are not equal. IEEE 754 fills the hole with subnormals — significands with a leading zero, allowed only at the lowest exponent. They give up significant digits to keep the spacing even all the way down to zero, and section 4 shows what that buys: two laws that hold with subnormals and fail without them.

Rounding is a function¶

The Handbook's second chapter frames rounding the way it is best understood. A rounding function ∘ sends every real number to a member, or to an infinity, and a floating-point operation is defined as a ⊕ b = ∘(a + b). The exact sum is taken first, as a real number, and only then rounded. IEEE 754 requires this correct rounding for +, −, ×, ÷ and √, which is why every float result on this page is the same on every machine that follows the standard.

There are five rounding functions, and section 3 prints all of them side by side:

RD, toward −∞: the largest member ≤ x.
RU, toward +∞: the smallest member ≥ x.
RZ, toward zero: RD for positive x, RU for negative x.
RN even, to nearest, with a tie going to the member whose significand is even. This is the default, and the only one Python's float ever uses.
RN away, to nearest, with a tie going to the larger magnitude.

The two RN functions agree except on an exact tie, and whatever RN returns is always RD(x) or RU(x); the program checks both on 1,025 inputs. Ties go to even by default because that has no bias: across many ties, about half go up and half go down, where ties-away pushes every one of them outward. Python's decimal module has all five under other names, and its ROUND_HALF_UP is ties-away, the rule most people learned at school.

Directed rounding is not a curiosity. RD and RU are how interval arithmetic ↗ keeps a lower and an upper bound that are guaranteed.

Which laws survive, and why those¶

Section 4 is a proof by exhaustion for the toy system: every pair and every triple of members was tried. In the output, (+) stands for ⊕. The pattern in the verdicts is the useful part.

law	verdict
a ⊕ b = b ⊕ a, and a ⊗ b = b ⊗ a	holds
a ≤ b implies a ⊕ c ≤ b ⊕ c	holds
a ⊖ b = 0 implies a = b	holds — only with subnormals
b/2 ≤ a ≤ 2b implies a ⊖ b is exact	holds — only with subnormals
(a ⊕ b) ⊕ c = a ⊕ (b ⊕ c)	fails
a ⊗ (b ⊕ c) = (a ⊗ b) ⊕ (a ⊗ c)	fails
b ≠ 0 implies a ⊕ b ≠ a	fails
a ⊗ (1 ⊘ a) = 1	fails

The laws that survive compare one exact result with itself. a + b and b + a are the same real number before any rounding happens, so they round to the same member. And a rounding function never reverses the order of two numbers, so monotonicity survives too.

The laws that fail compare two routes to a result, and each route rounds at different points along the way. (0.1 + 0.2) + 0.3 rounds the sum 0.30000000000000004 first and then rounds again. 0.1 + (0.2 + 0.3) adds two stored values whose errors happen to cancel to exactly 0.5, and rounds only at the last step. Neither answer is a mistake: each is the correctly rounded result of the operations actually performed. The sibling Rust library follows this to its consequence for a compiler in letting the compiler reorder a float sum ↗.

The fourth row is Sterbenz's lemma: when a and b are within a factor of 2 of each other, a − b is itself a member, so the subtraction involves no rounding at all. Section 5 confirms it on a real float — (0.1 + 0.2) - 0.3 returns exactly the difference of the two stored values. That is the precise form of what the next lesson, catastrophic cancellation, says in words: the subtraction does not create the error, it only exposes the error its inputs already carried.

Which numbers can be members at all¶

No amount of precision changes whether a fraction can be stored exactly. In radix β, a fraction in lowest terms has a finite expansion only when every prime factor of its denominator divides β. In radix 2 that means a power of two in the denominator, so 1/10 is out however many bits you allow. Radix 10 lets 1/10 in and still keeps 1/3 out, and radix 3 would do the reverse.

The Rust library walks through the binary expansion itself in what a float actually stores ↗. The point here is the one before it: the choice of β settles which fractions can ever be members, before a single digit is counted.

Off the ends of the range¶

A range limit is where a wrong answer can stay silent. Among its list of well-known arithmetic bugs, the Handbook records a release of Maple (7.0, from 2001) that returned 1 for 1001!/1000!, whose true value is 1001.

On IEEE 754 floats that computation cannot return 1. Both factorials overflow to inf long before 1000 — 171! already does — and inf / inf is nan, a value that announces itself. Section 6 shows the other possibility: a system that stopped at the largest float instead would compute largest / largest = 1.0, a plausible number that is completely wrong. That is the design reason IEEE 754 has infinities and NaNs at all.

The fixes are all ways of never leaving the set: cancel symbolically first, use exact integers, or work with logarithms, where 1001! becomes a number near 5,919 instead of a number with 2,571 digits.

The four numbers are the compromise¶

Section 7 puts four real formats side by side, and each parameter buys something and costs something:

Precision buys significant digits. Binary16 stores 1/10 wrong in the fourth significant digit, binary64 in the seventeenth. It costs memory and time.
Exponent range buys headroom. Binary16's largest value is 65,504, which is why 70,000 does not fit.
Radix 10 buys exact decimal fractions, which is why money uses it. It costs an uneven ruler — the relative gap varies tenfold across each power of 10, against twofold across each power of 2 — and, packed into 64 bits, fewer distinct values than binary64.
Radix 2 buys the most even ruler and the simplest hardware, and costs the decimal fractions people write every day.

Binary64 is the default not because it is correct — no finite set is — but because, for most numerical work, it is the best of these trades.

Run it yourself¶

From the root of your clone of this repository:

python3 01_Precision/machine_numbers/examples/machine_numbers.py

Auf Deutsch: Stichwörter¶

Eine Gleitkommazahl ist ein exaktes Element einer endlichen Menge; welche Rechengesetze das überlebt und welche nicht.

Stichwörter: Maschinenzahl, Gleitkommazahl (float, binary64), Mantisse, Exponent, Maschinengenauigkeit (machine epsilon), Rundungsfehler, Assoziativgesetz verletzt.