Floating-Point Arithmetic (IEEE 754)

A double is not a real number with a bit of fuzz on the end. It is an exact rational — a 53-bit integer times a power of two — drawn from a finite set of $1.844 \times 10^{19}$ values that are dense near zero and astronomically sparse near the top of the range. Every surprise in floating point follows from that one fact plus the rule that each operation returns the correctly rounded result of the exact operation on the exact inputs. IEEE 754 does not make arithmetic approximate; it makes the approximation reproducible, and this entry is about the machinery that gets you there.

The encoding #

A binary interchange format packs $k$ bits into three fields, in this order from the most significant bit:

 binary64 (double), k = 64
 ┌─┬───────────────┬────────────────────────────────────────────────────┐
 │S│  E (11 bits)  │            T (52 bits, trailing significand)        │
 └─┴───────────────┴────────────────────────────────────────────────────┘
  63 62           52 51                                                 0

The leading significand bit $d_0$ is not stored: it is implied by the exponent field, which is what buys the format one extra bit of precision. With $p$ the precision in bits, $t = p-1$ the stored significand width, $w$ the exponent width, and $E$ the biased exponent, IEEE 754-2019 §3.4 defines the value:

$$v = (-1)^S \times 2^{\,E-\text{bias}} \times \left(1 + 2^{\,1-p} \cdot T\right) \qquad \text{for } 1 \le E \le 2^w - 2$$

$$v = (-1)^S \times 2^{\,e_{\min}} \times \left(0 + 2^{\,1-p} \cdot T\right) \qquad \text{for } E = 0$$

Two exponent values are reserved and never mean "a normal number": $E = 0$ encodes zero and the subnormals, and $E = 2^w - 1$ encodes the infinities ($T = 0$) and the NaNs ($T \ne 0$). Everything in between is a normal number with an implicit leading 1.

Format parameters #

Values below are IEEE 754-2019 Table 3.5 (the parameters) with the derived limits computed from them. bfloat16 is not an IEEE interchange format — it is a vendor format that keeps binary32's exponent field and truncates the significand — but it is in the table because it is what most ML hardware actually stores.

Two rows deserve a second look. bfloat16 and binary32 have identical dynamic range — same $w$, same bias — which is the entire design argument: deep-learning training cares about range (gradients spanning many orders of magnitude) far more than about mantissa bits, so truncating binary32 to its top 16 bits is a free conversion (a shift) that never overflows where binary32 wouldn't. binary16 trades the other way: it has 3 more significand bits than bfloat16 but its max finite value is 65504, which real training runs hit.

And "double has 15–17 digits" is two different numbers. $53\log_{10}2 = 15.95$ is how many decimal digits are always reproduced correctly; 17 is how many you must print to guarantee that reading the text back recovers the identical bit pattern. Goldberg gives the same split for binary32: 9 digits round-trip a single. Printing 15 digits and parsing them back is a lossy operation.

Decoding a real bit pattern #

0.1 as a double is 0x3FB999999999999A:

 S = 0
 E = 0x3FB = 01111111011₂ = 1019  →  e = 1019 − 1023 = −4
 T = 0x999999999999A
 significand = 1.1001100110011001100110011001100110011001100110011010₂
 value       = (nearest double to 1.6) × 2⁻⁴
             = 0.1000000000000000055511151231257827021181583404541015625

The stored value is not 0.1. It cannot be: 1/10 in binary is the repeating fraction 0.0001100110011…, and any finite binary significand truncates it. The nearest double is above 0.1 by about $5.55 \times 10^{-18}$. This is the whole of the "floating point is inaccurate" folklore, and it is a base-conversion fact, not an arithmetic one — decimal and integer money types exist precisely because financial values are specified in base 10.

Normals, subnormals, and gradual underflow #

Without subnormals there would be a hole around zero $2^{52}$ times wider than the local spacing: the smallest normal double is $2^{-1022} \approx 2.225 \times 10^{-308}$, adjacent normals just above it are $2^{-1074}$ apart, and the next value below it would be zero. Subnormals fill that hole by giving up the implicit leading 1 and letting the significand slide right, trading precision for range one bit at a time down to $2^{-1074} \approx 4.94 \times 10^{-324}$.

Verify it directly — these are distinct subnormal doubles whose difference is representable:

double q = 0x1p-1070, r = 0x1p-1071;   // both subnormal
assert(q != r);
assert(q - r != 0.0);                   // 3.9525251667299724e-323, also subnormal

Precision: epsilon, ulp, and the sparsity of large numbers #

Watch the ulp grow. Bar length is logarithmic: $\text{width} = 8\log_{10}(\text{ulp}/2^{-52}) + 8$ pixels.

Above $2^{53} = 9\,007\,199\,254\,740\,992$ the ulp exceeds 1 and consecutive integers stop being representable: (double)(2^53 + 1) == (double)(2^53). This is the single most load-bearing number in practice.

Comparing with a fixed abs(a-b) < 1e-9 inherits this problem in the other direction: at magnitude $10^{12}$ that tolerance is smaller than one ulp, so it is an equality test; at magnitude $10^{-12}$ it accepts everything. Tolerances need to be relative, or expressed in ulps.

Rounding #

Every arithmetic operation in the standard is specified as: compute the exact result with unbounded range and precision, then round once to the destination format (§5.1). Five rounding-direction attributes exist (§4.3):

An implementation must provide roundTiesToEven plus the three directed attributes; roundTiesToAway is optional for binary. On x86 the mode lives in two bits of MXCSR (fesetround in C99); on AArch64 in FPCR.RMode.

Correct rounding is guaranteed for +, −, ×, ÷, sqrt, fma, remainder, and format conversions. It is not guaranteed for sin, exp, pow, or anything else in libm — accuracy there is implementation-defined, which is why the same program can produce different last bits on glibc versus musl versus the Apple libm.

Worked example: why 0.1 + 0.2 != 0.3 #

The usual explanation ("0.1 isn't exact") is only the first half. The interesting part is that the exact sum lands on a perfect tie, and ties-to-even pushes it away from 0.3.

Start with the two operands' significands, both 1.1001100110011001100110011001100110011001100110011010₂ (they are the same bits — 0.2 is 0.1 with the exponent one higher):

0.1 = 1.1001…1010₂ × 2⁻⁴      0x3FB999999999999A
0.2 = 1.1001…1010₂ × 2⁻³      0x3FC999999999999A

Align to $2^{-2}$ and add exactly. The exact sum is $10808639105689191 / 2^{55}$, which as a significand times $2^{-2}$ needs 54 bits:

exact sum = 1.0011001100110011001100110011001100110011001100110011|1₂ × 2⁻²
                                                     53 kept bits ^ ^ round bit = 1, sticky = 0

Round bit 1 with nothing after it is an exact halfway case. The two candidates:

…0110011₂ × 2⁻²  = 0.29999999999999998889776975374843…   0x3FD3333333333333   (LSB 1, odd)
…0110100₂ × 2⁻²  = 0.30000000000000004440892098500626…   0x3FD3333333333334   (LSB 0, even)

roundTiesToEven takes the even one. Meanwhile the literal 0.3 is a separate base-conversion, and its nearest double is the odd candidate, 0x3FD3333333333333. So:

Note what this rules out: it is not that the answer is "a bit off". It is off by the smallest possible amount, in a direction fixed by the tie-break rule, and it is bit-for-bit identical on every conforming machine on earth.

Signed zero, infinities, NaN #

Signed zero #

$+0$ and $-0$ are distinct bit patterns (0x0000000000000000 and 0x8000000000000000) that compare equal (§5.11: "Comparisons shall ignore the sign of zero"). The sign is not decoration — it carries the direction from which underflow arrived, which preserves identities that would otherwise break:

1.0 /  0.0  ==  +inf
1.0 / -0.0  ==  -inf         // so 1/(1/x) == x survives x = ±inf
sqrt(-0.0)  ==  -0.0
atan2( 0.0, -1.0) ==  +3.14159…
atan2(-0.0, -1.0) ==  -3.14159…   // branch cut of the complex log/sqrt

§6.3 also fixes the sign of an exact-zero sum: $x + (-x)$ is $+0$ under every rounding attribute except roundTowardNegative, where it is $-0$. And x + x and x - (-x) keep the sign of x under all of them.

Infinities #

Overflow under round-to-nearest produces $\pm\infty$ rather than the largest finite value: §4.3.1 specifies that any exact result of magnitude at least $2^{e_{\max}} \times (2 - 2^{-p})$ rounds to infinity. Infinity is a real value in the arithmetic, not an error state — 1/inf == 0, inf + 1 == inf, and inf - inf is NaN.

NaN #

A NaN is any $E = 2^w - 1$ pattern with $T \ne 0$. §6.2.1: the leading bit of $T$ (that is, $d_1$) distinguishes the two kinds — 1 for quiet, 0 for signalling — and the remaining $t-1$ bits are the payload, free for diagnostic use. For binary64 that means $2^{51}$ quiet and $2^{51}-1$ signalling patterns per sign bit: roughly $9 \times 10^{15}$ distinct NaNs, of which one, 0x7FF8000000000000, is the canonical quiet NaN most hardware produces.

Signalling NaNs exist for uninitialised-variable detection: fill memory with sNaNs and the first arithmetic use traps. In practice almost nothing uses them, because most languages quiet an sNaN as soon as a value is copied through a register on some platforms.

Non-associativity #

Floating-point addition is commutative but not associative, because each + rounds:

(1e16 + -1e16) + 1.0  ==  1.0
 1e16 + (-1e16 + 1.0) ==  0.0     // the 1.0 is absorbed: ulp(1e16) = 2

(1.0 + 1e100) + -1e100 ==  0.0    // the 1.0 is absorbed
 1.0 + (1e100 + -1e100) ==  1.0

Nor is it distributive, and x/10.0 is not x*0.1 (the second rounds twice: once converting the literal, once multiplying). The direct consequences:

• Parallel reductions are not reproducible. A summation split across 8 threads and one split across 16 give different last bits, and re-running with a different thread count changes the answer. If bit-reproducibility matters, you need a fixed reduction tree, a deterministic order, or a compensated/exact accumulator — see rounding error and cancellation.
• Vectorisation is a semantic change. Turning a serial sum into four SIMD lanes reassociates it. Compilers will not do this under plain -O2/-O3 for exactly that reason; -ffast-math / -Ofast / #pragma omp simd reduction explicitly grant permission.
• Summation order matters. Adding a large accumulator to many small terms loses the small ones one at a time; sorting ascending, pairwise summation, or Kahan compensation all fix it to different degrees.

FMA: one rounding instead of two #

fusedMultiplyAdd(x, y, z) computes $(x \times y) + z$ "as if with unbounded range and precision, rounding only once to the destination format" (§5.4.1). Two consequences follow, and they pull in opposite directions.

It is more accurate. The full $2p$-bit product survives into the addition instead of being truncated to $p$ bits first. A discriminant $b^2 - 4ac$ with $a = c = 1$, $b = 2 + 2^{-31}$:

double b = 2.0 + 0x1p-31, w = 4.0;
b*b - w        // 1.8626451492309570312e-09  == 2⁻²⁹ exactly — the 2⁻⁶² tail was lost in b*b
fma(b, b, -w)  // 1.8626451494477974657e-09  == 2⁻²⁹ + 2⁻⁶² — the exact answer

The naive form is wrong by a relative $2^{-33} \approx 1.2 \times 10^{-10}$ — about six orders of magnitude worse than $\varepsilon$ — because $b^2$ rounded away exactly the bits that the subtraction was going to expose. This is why FMA is the standard tool for determinants, discriminants, dot products, polynomial evaluation (Horner with FMA is both faster and more accurate), and Newton iterations for div/sqrt.

It gives you the error term for free. Since the product's rounding error is itself representable, one FMA extracts it exactly:

double p = x*y;
double e = fma(x, y, -p);   // p + e == x*y exactly, in double-double form

That is TwoProduct, the multiplication half of every double-double / compensated-arithmetic library. Verified: fma(0.1, 3.0, -(0.1*3.0)) is exactly $-2^{-55}$.

Subnormals cost real cycles #

Subnormals are handled by a slower path — sometimes a microcode assist — on essentially all x86 silicon. Wittmann et al. measured AVX double-precision throughput in cycles per scalar operation with FTZ+DAZ enabled ("F+D") versus disabled:

Read the numbers, not the vibe. Overflow and NaN operands cost nothing on these cores — a NaN-poisoned array runs at full speed. The penalty is specifically subnormal inputs and underflowing results, and it is roughly two orders of magnitude: 0.25 → ~32–40 cycles for add/mul, 6.7 → ~57 for divide. Intel improved it across generations (38.2 → 37.7 → 31.9 cycles for an underflowing add from Sandy Bridge to Haswell), but did not remove it. Curiously, an underflow in the multiply half of an FMA costs nothing on Haswell, while an underflow in the addition half costs 33.5 cycles.

Two mode bits in MXCSR opt out: FTZ (bit 15) returns zero instead of a subnormal result, DAZ (bit 6) treats subnormal inputs as zero. Neither is IEEE 754 conforming. AArch64 has the equivalent in FPCR.FZ.

Practical rules #