Floating-point arithmetic and precision

Why it exists

The page on IEEE 754 tells how a number is stored. Doing arithmetic with the stored numbers is a different matter: the sum of two floats is usually not a float. $2^{24} + 1$ needs 25 significant bits and a single-precision number has 24. The hardware must therefore produce a nearby float, and the rules for which one are fixed by the standard so that every machine gives the same answer. Exam questions ask for the bit pattern of a sum, the number of bits needed, or why a program prints a surprising result; the official syllabus asks for "number representation and arithmetic (fixed and floating point)".

Prerequisites

You need the fields of a single-precision number, the hidden bit, the bias, the special values and denormals from IEEE 754, and binary addition from adders and subtractors.

Precise definitions

The rule of the standard. For each of $+, -, \times, \div$ and the square root, the computed result must be the exact mathematical result rounded to a float. Rounding means choosing the nearest float; if the exact result is exactly halfway between two floats, the one whose last fraction bit is 0 (the even one) is chosen: round to nearest, ties to even, the default mode. The other IEEE modes (toward zero, toward $+\infty$, toward $-\infty$) are not computed here. A result too large for the format becomes $\pm\infty$; one below the smallest normal number becomes a denormal (gradual underflow) or zero.

Spacing. In $[2^k, 2^{k+1})$ the floats are $2^{k-23}$ apart in single precision (one ulp, unit in the last place). The relative gap is between $2^{-24}$ and $2^{-23}$. $2^{-23}$ is the machine epsilon (the gap just above 1). The relative error of one rounded operation (normal results) is at most half a gap divided by the number, $2^{-24}$.

Guard, round and sticky bits. To obtain the correctly rounded result the adder keeps three extra bits to the right of the significand: G (guard, the first bit beyond), R (round, the next), and S (sticky, the OR of all bits beyond those). Three bits are enough to round correctly; the sticky bit is what makes "exactly halfway" different from "just above halfway".

Special values in arithmetic. Any operation with a NaN gives NaN. $\infty + \infty = \infty$, $\infty + (-\infty) =$ NaN, $\infty \times 0 =$ NaN, $x + \infty = \infty$, $x \times \infty = \pm\infty$ for $x \ne 0$. $x + (-x) = +0$ (in the default mode). A comparison with NaN is false, even $\text{NaN} = \text{NaN}$.

How it works

Addition and subtraction

1. Unpack both numbers: sign, exponent, significand with the hidden bit.
2. Compare exponents. The number with the smaller exponent is the one that is shifted. The difference $d$ is the shift.
3. Align: shift the smaller significand right by $d$ places, keeping G, R and S for what is shifted out.
4. Add the significands if the signs are equal, subtract if they differ.
5. Normalise: a carry out needs a shift right by 1 (exponent $+1$); cancelled leading bits need a shift left (exponent decreased by the same amount).
6. Round to nearest even using the last kept bit, G, R and S.
7. Pack (and handle overflow, a rounding carry that needs another normalisation, and underflow).

$12.5 + 0.1875$: the exponents differ by 6, the smaller number's significand slides six places to the right, the sum fits in the 24 bits, nothing needs to be rounded. The result 12.6875 is exact.

Rounding to nearest even

If the guard bit is 0 the discarded part is below half a place: truncate. If G is 1 and anything else beyond it is 1, it is above half: round up. Only G=1 with R=S=0 is a tie; then the kept bits are rounded to the even neighbour, up if the last kept bit is 1, unchanged if it is 0. Over many operations this does not drift in one direction, which is why it is the default and why it is not "round half up".

$2^{24} + 3$: near $2^{24}$ the floats are 2 apart, so $2^{24} + 3$ lies exactly halfway between $2^{24} + 2$ and $2^{24} + 4$. The significand after the alignment is $\dots 0001 | 100$: G=1, R=S=0, last kept bit 1 (odd): round up to the even neighbour, $2^{24} + 4$. By the same rule $2^{24} + 1$ (halfway between $2^{24}$ and $2^{24} + 2$, last kept bit 0) rounds down to $2^{24}$: adding 1 to $2^{24}$ changes nothing, and $2^{24}$ is the first integer whose successor is not a single-precision float.

Cancellation

When two nearly equal numbers are subtracted, the leading bits cancel; the result is computed exactly (nothing is rounded: the exponent difference is 0) but the normalisation shifts left by 23 places and the result has only one significant bit. If $x$ and $y$ themselves carried rounding errors in their last bits, those bits are now the whole result: the relative error is huge. This is catastrophic cancellation, the reason numerical code avoids subtracting nearly equal numbers.

Multiplication

Multiplication needs no alignment. Exponents add (and one bias is subtracted, because both fields contain it), significands multiply, the signs xor. The product of two numbers in $[1, 2)$ is in $[1, 4)$: if it is at least 2, shift right by one and add one to the exponent. Then round.

How big are the gaps?

The gap is $2^{k-23}$ for a number in $[2^k, 2^{k+1})$: it doubles every time the exponent grows by one. A number smaller than half the gap, added to $x$, is rounded away completely: this is absorption. All integers up to $2^{24} = 16\,777\,216$ are single-precision numbers (the gap is 1 or less); above it, only even numbers (up to $2^{25}$), then multiples of 4, and so on. There are $2^{23} = 8\,388\,608$ floats in $[1, 2)$, and the same number in every other interval $[2^k, 2^{k+1})$.

Addition is not associative

Each addition is rounded, so what is lost depends on the order. $(2^{24} + 1) + 1 = 2^{24}$ (each 1 is absorbed), but $2^{24} + (1 + 1) = 2^{24} + 2$. In the second row each $2^{-24}$ is exactly half a gap below 1, a tie that rounds back to 1 twice, while $2^{-24} + 2^{-24} = 2^{-23}$ is a whole gap and survives. In the last row $(10^8 - 10^8) + 1 = 1$ but $10^8 + (-10^8 + 1) = 0$, because $-10^8 + 1$ rounds to $-10^8$. Even for random numbers in $[0, 1)$ about a quarter of all triples give different results in the two orders. Addition and multiplication are commutative; they are not associative or distributive, and so a compiler may not reorder floating-point sums unless it is told that exactness does not matter.

A real program

The program is compiled with plain float arithmetic and prints what the hardware does. volatile stops the compiler from computing the sums at compile time.

#include <float.h>
#include <stdint.h>
#include <stdio.h>
#include <string.h>
 
static uint32_t bits(float x) {
    uint32_t u;
    memcpy(&u, &x, sizeof u);
    return u;
}
 
int main(void) {
    volatile float big = 16777216.0f;          /* 2^24 */
    volatile float one = 1.0f, three = 3.0f, tiny = 5.9604645e-8f;   /* 2^-24 */
    float s1 = big + one, s3 = big + three;
    printf("2^24 + 1 = %.1f   (bits %08X)\n", s1, bits(s1));
    printf("2^24 + 3 = %.1f   (bits %08X)\n", s3, bits(s3));
    float left = (one + tiny) + tiny, right = one + (tiny + tiny);
    printf("(1 + 2^-24) + 2^-24 = %.9g   1 + (2^-24 + 2^-24) = %.9g\n", left, right);
    volatile float a = 100000000.0f, b = -100000000.0f;
    float l2 = (a + b) + one, r2 = a + (b + one);
    printf("(1e8 - 1e8) + 1 = %.1f   1e8 + (-1e8 + 1) = %.1f\n", l2, r2);
    printf("FLT_EPSILON = %.9g (this is 2^-23: %s)\n", FLT_EPSILON, FLT_EPSILON == 1.1920928955078125e-07f ? "yes" : "no");
    volatile double d1 = 0.1, d2 = 0.2, d3 = 0.3;
    printf("double: 0.1 + 0.2 == 0.3 is %s, 0.1 + 0.2 = %.17g\n", d1 + d2 == d3 ? "true" : "false", d1 + d2);
    return 0;
}

2^24 + 1 = 16777216.0   (bits 4B800000)
2^24 + 3 = 16777220.0   (bits 4B800002)
(1 + 2^-24) + 2^-24 = 1   1 + (2^-24 + 2^-24) = 1.00000012
(1e8 - 1e8) + 1 = 1.0   1e8 + (-1e8 + 1) = 0.0
FLT_EPSILON = 1.1920929e-07 (this is 2^-23: yes)
double: 0.1 + 0.2 == 0.3 is false, 0.1 + 0.2 = 0.30000000000000004

Worked examples

Cheat box

Traps and GATE patterns

Connections

The adder of the significands is the binary adder of the digital design pages, with a shifter in front and a normaliser after it. Fixed-point arithmetic has no alignment step but a fixed range; floating point trades that for a range of about $10^{\pm 38}$ and a relative accuracy of $10^{-7}$. The C types float and double of C types and operators use these formats.

PYQs

Consider the real valued variables X, Y and Z
represented using the IEEE -754 single precision
floating-point format. The binary representations
of X and Y in hexadecimal notation are as
follows:

X: 35C00000
Y: 34A00000

Let Z = X + Y.

Which one of the following is the binary
representation of Z, in hexadecimal notation?

Self-check

Revision card

References