🎯Computer Arithmetic and Error

Floating-Point Representation

Binary fractions, IEEE 754, machine epsilon, rounding, overflow and underflow

Dhahran, 25 February 1991

Saudi Arabia, the last week of the Gulf War. A Patriot air-defence battery at Dhahran has been switched on for about 100 hours. The Patriot predicts where an incoming missile will be next by multiplying its speed by elapsed time. Its computer counts time in tenths of a second, as an integer, and converts to seconds by multiplying by 0.1 held in a 24-bit register. But 0.1 has no finite binary expansion. In binary it is 0.0001100110011001100… forever. The register keeps 23 bits after the point and throws the rest away. Each tick is short by about 0.000000095 seconds. Nobody would notice that in an hour. After 100 hours — 3.6 million ticks — the clock is 0.34 seconds behind. A Scud travels well over half a kilometre in 0.34 s. When a Scud came in that evening, the radar looked for it where the drifted clock said it should be, found nothing in its range gate, and did not engage. The missile hit a barracks. 28 American soldiers were killed. The fault was known: a software fix that corrected the clock reached Dhahran the next day. This lesson is about what a computer actually stores when you write 0.1, and why the answer is never quite 0.1.

A real number can carry infinitely many digits. A computer word has 32 or 64 bits. Every number a program uses is therefore an approximation chosen from a finite set, and the whole of numerical analysis starts from one question: how far is the stored number from the one you meant, and what happens to that gap as the computation proceeds?

Integers are easy — every integer up to

2^{53}

fits exactly in a double. Fractions are not. In base 10 the fraction 1/3 never terminates; in base 2 the same is true of every fraction whose denominator is not a power of two, including 1/10. So 0.1, 0.2 and 0.3 are all approximations in binary, and

0.1 + 0.2

prints as

0.30000000000000004

Derivation

Why 0.1 never ends in binary

To write a fraction in base 2, multiply it by 2 again and again. The whole-number part that appears each time is the next binary digit; the fractional part is carried on to the next line.

0.1 \times 2 = 0.2 \;\to\; 0, \qquad 0.2 \times 2 = 0.4 \;\to\; 0, \qquad 0.4 \times 2 = 0.8 \;\to\; 0,

0.8 \times 2 = 1.6 \;\to\; 1, \qquad 0.6 \times 2 = 1.2 \;\to\; 1, \qquad 0.2 \times 2 = 0.4 \;\to\; 0, \;\dots

The fractional part

0.2

has come back, so from here the same four digits

0011

repeat for ever:

0.1 = 0.0\,0011\,0011\,0011\ldots_2

Why it must repeat: every fractional part that appears has the form

k/10

with

k

between 0 and 9, so after at most ten doublings one of them returns. The expansion stops only if a fractional part becomes exactly 0, which needs

2^n \cdot \tfrac{1}{10}

to be a whole number for some

n

. Because

10 = 2 \times 5

, no power of 2 cancels the 5. The same argument shows that a fraction terminates in binary exactly when its denominator is a power of two.

Scientific notation fixes the precision, not the size:

6.022 \times 10^{23}

and

1.6 \times 10^{-19}

both have four significant digits. Floating point is scientific notation in base 2. A number is stored as a sign, an exponent and a fixed number of significant bits, so the relative precision is the same whether the number is the mass of a galaxy or the charge of an electron.

x = (-1)^{s}\,\left(1.f_1 f_2 \dots f_{52}\right)_2 \times 2^{\,e-1023}

IEEE 754 binary64 (double). s is one bit, the biased exponent e is 11 bits, the fraction f is 52 bits. Because a normalised binary number always starts with 1, that leading 1 is not stored — 53 significant bits for the price of 52.

Format	Sign / exponent / fraction bits	Machine epsilon	Largest	Smallest normal	Decimal digits
binary32 (single)	1 / 8 / 23	$2^{-23} \approx 1.19 \times 10^{-7}$	$3.40 \times 10^{38}$	$1.18 \times 10^{-38}$	about 7
binary64 (double)	1 / 11 / 52	$2^{-52} \approx 2.22 \times 10^{-16}$	$1.80 \times 10^{308}$	$2.23 \times 10^{-308}$	about 16

The two IEEE 754 formats every engineering code uses

💡Machine epsilon and unit roundoff

Machine epsilon

\varepsilon_{mach}

is the gap between 1 and the next larger float:

2^{-52}

in double precision. With round-to-nearest, any real

x

in range is stored as

\mathrm{fl}(x) = x(1+\delta)

with

|\delta| \le u

, where the unit roundoff is

u = \varepsilon_{mach}/2 = 2^{-53} \approx 1.1 \times 10^{-16}

. This one inequality is the model of every error analysis in this course: storing a number costs a relative error of at most

u

. Not an absolute error — the gaps between floats grow with the numbers. Between 1 and 2 they are

2.2 \times 10^{-16}

apart; between

2^{52}

and

2^{53}

they are 1 apart; above

2^{53}

the odd integers are not representable at all.

Derivation

Why rounding costs at most

u = \varepsilon_{mach}/2

Take any real

x

between two powers of two,

2^e \le x < 2^{e+1}

. Every double in that range has the same exponent

e

and differs only in its 52 fraction bits, so neighbouring doubles are a fixed distance apart: the value of the last bit,

\text{gap} = 2^{-52} \times 2^{e} = \varepsilon_{mach}\, 2^{e}.

Round-to-nearest moves

x

to the closer of its two neighbours, so it moves it by at most half a gap:

|\mathrm{fl}(x) - x| \le \tfrac12\,\varepsilon_{mach}\,2^{e}.

Divide by

|x|

, which is at least

2^e

\frac{|\mathrm{fl}(x) - x|}{|x|} \le \frac{\tfrac12\varepsilon_{mach}2^{e}}{2^{e}} = \frac{\varepsilon_{mach}}{2} = u.

Call that relative error

\delta = (\mathrm{fl}(x) - x)/x

; then

\mathrm{fl}(x) = x(1 + \delta)

with

|\delta| \le u

— the model above. Two consequences fall out of the same lines. The absolute gap

\varepsilon_{mach}2^e

doubles with every power of two, which is why it is 1 between

2^{52}

and

2^{53}

. And chopping can move

x

by up to a whole gap, always towards zero, so its bound is

\varepsilon_{mach}

, not

\varepsilon_{mach}/2

📐What is stored for 0.1

Double precision rounds 0.1 to the nearest float, which is exactly

0.1000000000000000055511151231257827021181583404541015625

The error is

5.55 \times 10^{-18}

, a relative error of

5.55 \times 10^{-17}

— below

u = 1.1 \times 10^{-16}

, exactly as the model promises. In the binary pattern the exponent is

-4

(because

0.1 = 1.6 \times 2^{-4}

) and the fraction is the repeating pattern 1001 1001 … rounded up in its last bit. The Patriot's register did two things worse: it kept only 23 fraction bits, and it chopped instead of rounding. Chopping 0.1 at 23 bits gives 0.09999990463… — an error of

9.54 \times 10^{-8}

, about a billion times larger.

Derivation

Where the Patriot's

9.54 \times 10^{-8}

s comes from

From the doubling table above, the binary digits of 0.1 are

b_1 b_2 b_3 \ldots = 0\,0011\,0011\,0011\ldots

: digit

b_k

is 1 when

k - 2

leaves remainder 2 or 3 on division by 4.

The register kept

b_1

b_{23}

and threw the rest away. The discarded tail starts at

b_{24}

, and

24 - 2 = 22

leaves remainder 2, so the tail reads

1100\,1100\ldots

0.1 - \text{chopped} = \sum_{k \ge 24} b_k 2^{-k} = 2^{-23} \times (0.1100\,1100\ldots)_2.

A block of four digits repeated for ever is worth the block divided by

2^4 - 1 = 15

(the binary version of

0.\overline{12} = 12/99

in decimal). The block

1100_2

is 12, so

(0.1100\,1100\ldots)_2 = 12/15 = 0.8

. Therefore

0.1 - \text{chopped} = 0.8 \times 2^{-23} = 9.537 \times 10^{-8}.

The clock counted ticks correctly; each tick was converted to seconds with a number that was too small by this amount, always in the same direction.

The purpose of computing is insight, not numbers.

– Richard Hamming, Numerical Methods for Scientists and Engineers (1962)

Rounding modes. IEEE 754 rounds to nearest by default, breaking ties to the value whose last bit is even. The alternatives — toward zero (chopping), toward

+\infty

, toward

-\infty

— exist for interval arithmetic and special uses. Round-to-nearest errs in both directions, so independent errors partly cancel. Chopping always errs the same way, so its errors add up linearly. That difference is the whole Patriot story:

9.54 \times 10^{-8}\ \text{s per tick} \times 3.6 \times 10^{6}\ \text{ticks} = 0.343\ \text{s}

The US General Accounting Office's report tabulated exactly this growth: 0.0275 s after 8 hours, 0.1648 s after 48 hours, 0.3433 s after 100 hours, by which point the range gate had shifted by 687 m.

Running out of range. The exponent has limits too.

Overflow — a result larger than about $1.8 \times 10^{308}$ becomes $\infty$ . Squaring $10^{200}$ overflows even though the answer you want, say $\sqrt{a^2+b^2}$ , is modest; math.hypot exists precisely to avoid that.
Underflow — results smaller than $2.2 \times 10^{-308}$ enter the subnormal range, where the leading 1 is given up so the gap to zero is filled evenly. Precision is lost gradually; below $4.9 \times 10^{-324}$ the result is 0.
Inf and NaN — $1/0 = \infty$ , $\infty - \infty =$ NaN (not a number). NaN compares unequal to everything, itself included, and silently poisons every calculation it touches.
Integer conversion — casting a float into a fixed-size integer can overflow too, and that is not always caught.

🛰Ariane 5, Flight 501 — 4 June 1996

Thirty-seven seconds after lift-off from Kourou, Ariane 5's first flight veered off course and broke up. The inertial reference software, reused from Ariane 4, converted a 64-bit floating-point value (the horizontal bias) into a 16-bit signed integer. Ariane 5 accelerated harder than Ariane 4, the value exceeded 32,767, and the conversion raised an operand error. The protection for that conversion had been removed to save processor time, on the grounds that Ariane 4's trajectory could never trigger it. Both inertial units failed the same way, within 72 milliseconds of each other. The inquiry board's report is short and worth reading in full.

💡Where you will use this

Every later course runs on doubles. GNSS Navigation: a pseudorange of 20,000 km stored with relative error

10^{-16}

is still good to 2 nm — but a receiver clock in single precision would not be. Orbital Mechanics and Attitude Control: flight software often runs in single precision to save power, so its developers do this lesson's error budget on every variable. Rocket Propulsion and Structures: nothing breaks until a value overflows or a NaN propagates through a solver at 3 a.m.

📚Further reading

David Goldberg, What Every Computer Scientist Should Know About Floating-Point Arithmetic (1991) — the classic tutorial, free online.
Michael Overton, Numerical Computing with IEEE Floating Point Arithmetic (SIAM, 2001) — a short book on exactly this lesson.
Nicholas Higham, Accuracy and Stability of Numerical Algorithms (2nd ed., SIAM, 2002), chapter 2, 'Floating Point Arithmetic'.
US General Accounting Office, Patriot Missile Defense: Software Problem Led to System Failure at Dhahran (IMTEC-92-26, 1992).

↓But there is a problem

So every stored number carries a relative error of at most about

10^{-16}

. That sounds like the end of the story: sixteen digits is far more than any measurement deserves. But programs do not store one number. They add, subtract, multiply and divide millions of them, and each operation is rounded again. In 1982 a stock exchange in Vancouver launched a new index, and every operation on it lost less than a thousandth of a point. Twenty-two months later the index had lost half its value — while the market it measured had gone up.