Chapter 20: Floating-Point Arithmetic — Precision, Rounding, and Why 0.1 + 0.2 ≠ 0.3

"The first principle is that you must not fool yourself — and you are the easiest person to fool." — Richard Feynman, "Cargo Cult Science" (1974)

Overview

Open a Python prompt, or any language with an interactive console, and type 0.1 + 0.2. The answer is not 0.3. It is 0.30000000000000004. This is not a bug in Python, and it will not be a bug in your Fortran either — it is the single most important fact about computer arithmetic, and every scientific result you ever compute rests on making peace with it. The real numbers you learned about in school are infinite and continuous; the numbers your processor actually stores are finite and discrete, a sparse grid of representable values with gaps between them. Most of the decimals you write down — including the utterly ordinary 0.1 — fall into a gap and are silently replaced by the nearest number on the grid. The arithmetic that follows is the arithmetic of the grid, not of the real line, and the difference, usually invisible, occasionally destroys an answer completely.

This chapter opens Part V, the numerical heart of the book, and it comes first for a reason: linear algebra, integration, differential equations, and the heat solver you have been building all rest on floating-point arithmetic, and a method that ignores its properties will produce numbers that look plausible and are wrong. By the end of the chapter you will understand exactly what a real(dp) is down to its bits, you will be able to predict and bound the rounding error in a calculation, you will recognize the one arithmetic operation — subtraction of nearly equal numbers — that turns tiny errors into catastrophic ones, and you will know how to detect the NaN and Inf values that signal a computation has gone off the rails. This is the chapter that turns you from someone who uses floating-point numbers into someone who understands them.

In this chapter, you will learn to:

• Read the anatomy of an IEEE 754 real — sign, exponent, and significand — and say precisely what single, double, and quadruple precision can and cannot represent.
• Put numbers on the precision of a kind using epsilon, spacing, nearest, huge, and tiny, and explain why the gap between representable numbers grows with their magnitude.
• Spot catastrophic cancellation and rearrange the algebra so it never happens.
• Reason about overflow, underflow, NaN, and Inf, and detect them with the ieee_arithmetic module.
• Tell a well-conditioned problem from an ill-conditioned one, and a stable algorithm from an unstable one — the intuition that governs every method in the rest of Part V.
• Choose sp, dp, or qp on purpose, trading accuracy against memory and speed like the computational scientist you are becoming.

Learning Paths

How to read this chapter by track. - 🔬 Scientist ("I need my numbers to be right") — this is one of the most important chapters in the book for you. Read §20.2 (how big is the error), §20.3 (cancellation, the error you cause yourself), and §20.5 (conditioning and stability) closely. They are the difference between a trustworthy result and a plausible lie. - 📖 Standard — read straight through; every later numerical chapter assumes this material. - 🔧 Legacy ("I inherited old code") — §20.1 explains the double precision you will see everywhere, and §20.6 explains the accuracy consequences of an old code's single-precision choices. §20.4's ieee_arithmetic is a modern tool for hunting the NaNs that legacy codes so often emit. - ⚡ HPC ("I need it fast") — §20.6 is your chapter: single precision is often twice as fast for a bandwidth-bound kernel, and knowing when you can afford it is a real performance lever. The ⚡ Performance Note in §20.2 and the Project Checkpoint quantify the trade.

20.1 IEEE 754: How a Real Is Stored

In Chapter 3 we defined dp as selected_real_kind(15, 307) and made a rule of it — every real in this book is real(dp), every literal carries a _dp suffix — but we deliberately deferred the question of why a real is only approximate, and what those numbers 15 and 307 really mean. This is the chapter that pays that debt. A real(dp) value is not an abstract number; it is 64 specific bits arranged in a specific way, and once you can read those bits, the whole subject stops being mysterious.

Almost every computer built in the last thirty-five years represents real numbers according to a single standard, and it is worth naming it precisely because its guarantees are what let your Fortran produce the same answer on your laptop, a colleague's Mac, and a supercomputer in another country.

Definition (IEEE 754). IEEE 754 is the international standard for binary floating-point arithmetic, first published in 1985 (revised in 2008 and 2019). It specifies how real numbers are encoded in bits, how many bits each format uses, how the basic operations round their results, and what happens in exceptional cases like division by zero. Its central achievement is reproducibility: because the standard pins down the encoding and the rounding, a conforming calculation gives bit-for-bit identical results across conforming machines. gfortran's real kinds are IEEE 754 formats, which is why the selected_real_kind(15, 307) you wrote in Chapter 3 lands exactly on the standard's 64-bit format on every mainstream system.

A floating-point number is scientific notation in binary. Just as you can write a decimal number as $\pm d.dddd \times 10^{e}$ — a sign, a fractional part called the significand (or mantissa), and an exponent — IEEE 754 writes a number in base 2 as

$$ x = (-1)^{s} \times 1.f \times 2^{e} $$

and stores the three pieces in three adjacent bit-fields. For double precision (the 64-bit format, Fortran's real64 and our dp) the split is:

 double precision (binary64), 64 bits total:

 ┌─┬───────────────┬──────────────────────────────────────────────────┐
 │s│  exponent (11)│                 fraction (52)                     │
 └─┴───────────────┴──────────────────────────────────────────────────┘
  1       11                              52

   sign     biased exponent e+1023        the bits after the "1." of 1.f

• 1 sign bit s: 0 for positive, 1 for negative.
• 11 exponent bits: the exponent $e$, stored with a bias of 1023 (the stored value is $e + 1023$), so $e$ ranges from $-1022$ to $+1023$ for ordinary numbers.
• 52 fraction bits: the bits of the significand after the binary point. Here is the clever part: because a normalized binary number always starts with a 1 (there is no other nonzero digit in base 2), that leading 1 is not stored — it is implicit. So 52 stored bits buy you 53 bits of significand precision. This is the famous "hidden bit."

Definition (significand). The significand (also called the mantissa) is the part of a floating-point number that carries its significant digits — the $1.f$ in $\pm 1.f \times 2^{e}$. Its width sets the number's precision: double precision's 53-bit significand corresponds to about $53 \times \log_{10}2 \approx 15.95$ decimal digits, which is why we say double gives roughly 15–16 significant decimal figures. The exponent, by contrast, sets the number's range, not its precision.

Those two roles — significand for precision, exponent for range — are the key to reading the arguments of selected_real_kind(15, 307). The 15 asked for 15 decimal digits of precision, which needs the 53-bit significand; the 307 asked for a decimal exponent range out to $10^{307}$, which needs the 11-bit exponent. The standard offers three sizes, and you should know all three:

Now we can finally see why 0.1 is not 0.1. In base 2, the decimal fraction $0.1$ is a repeating fraction, exactly the way $1/3 = 0.3333\ldots$ repeats in base 10:

$$ 0.1_{10} = 0.0001100110011001100\ldots_2 \quad (\text{the } 1100 \text{ repeats forever}) $$

A 53-bit significand cannot hold infinitely many bits, so the value is chopped and rounded to the nearest representable number, which is not exactly $0.1$ but

$$ 0.1000000000000000055511151231257827\ldots $$

— a hair too large. The same happens to 0.2, also a hair too large. Add the two rounded values and the two little excesses add too, pushing the sum just past the rounded value of 0.3, which happens to sit a hair the other side. The result is a number that is not the stored 0.3, and Fortran, asked whether they are equal, correctly answers no.

program point_one_plus_point_two
  use kinds, only: dp
  implicit none
  real(dp) :: a = 0.1_dp, b = 0.2_dp, c = 0.3_dp

  print '(a, f20.17)', '0.1 + 0.2            = ', a + b
  print '(a, f20.17)', '0.3                  = ', c
  print '(a, l1)',     '(0.1 + 0.2) == 0.3 ? ', (a + b == c)
  print '(a, es13.6)', '(0.1 + 0.2) - 0.3    = ', (a + b) - c
end program point_one_plus_point_two

$ gfortran -std=f2018 -Wall point_one_plus_point_two.f90 -o ppp && ./ppp
0.1 + 0.2            =  0.30000000000000004
0.3                  =  0.29999999999999999
(0.1 + 0.2) == 0.3 ? F
(0.1 + 0.2) - 0.3    =  5.551115E-17

Read every line of that output slowly, because it overturns a lifetime of intuition. The sum 0.1 + 0.2 displays as 0.30000000000000004 — the true stored value is $0.3000000000000000444\ldots$, and f20.17 rounds it to seventeen decimals for display. The literal 0.3 displays as 0.29999999999999999 — it, too, is not exactly three-tenths, merely the closest double. The two differ by 5.551115E-17, which is not noise: it is exactly $2^{-54}$, the distance between two adjacent doubles near $0.3$. (One processor-dependent detail: whether a 0 is printed before the decimal point under the F descriptor is left to the compiler by the standard; gfortran prints it, so you see 0.30000...; some compilers print .30000....)

🚪 Threshold Concept — floating-point numbers are not the real numbers. This is the idea that changes everything after it. A real(dp) is not a real number; it is one of about $2^{64}$ specific rational numbers, spaced unevenly along the number line, with nothing at all in between. Every literal you write, every result you compute, is snapped to the nearest one. Arithmetic on these values obeys its own laws, not the laws of the real numbers: it is commutative but not associative, most decimals are not representable, and equality is a trap. Once you stop picturing a smooth continuum and start picturing a finite, gap-riddled grid, the rest of numerical computing — why we bound errors instead of expecting exactness, why we never test floats with ==, why the order of a summation matters — follows naturally. Hold this picture; the whole of Part V leans on it.

⚠️ Common Pitfall — never compare floats with ==. Because 0.1 + 0.2 is not the stored 0.3, writing if (x == 0.3_dp) is almost always a bug: it asks whether two computed values landed on the exact same grid point, which rounding rarely arranges. Test with a tolerance instead — if (abs(x - 0.3_dp) < 1.0e-12_dp) — or, better, a tolerance scaled to the magnitudes involved, as §20.2 will make precise. The one time == is legitimate is comparing against a value you know is exactly representable, such as 0.0_dp used purely as a sentinel.

🐍 Python Comparison. If you have met 0.1 + 0.2 == 0.30000000000000004 in Python, you have already met this chapter — because Python's float is exactly IEEE 754 double precision, the very same 64-bit binary64 as Fortran's real(dp). The bits are identical; only the printing differs. Python's repr shows the shortest decimal that round-trips (0.30000000000000004), while Fortran shows whatever your format descriptor asks for, but underneath they are the same number, computed by the same hardware. Two genuine differences are worth noting. First, Python's integers are arbitrary-precision — 10**400 is exact in Python and overflows a Fortran integer — but Python's floats are not, and inherit every quirk in this chapter. Second, NumPy's float64 is the same double, so a NumPy array of float64 and a Fortran real(dp) array hold bit-identical values, which is exactly why the two interoperate so cleanly through f2py (Chapter 15). The lesson: floating-point behavior is a property of the hardware standard, not of the language, so everything you learn here transfers everywhere.

📜 From History. Before IEEE 754, every vendor had its own floating-point format, with its own rounding, its own range, and its own idea of what to do on overflow. The same program gave different answers on an IBM mainframe, a DEC VAX, and a Cray, and porting numerical code was a nightmare of silent inaccuracies. William Kahan led the design of the 1985 standard — work for which he received the Turing Award — and its adoption is one of the quiet triumphs of computing: the reason your simulation reproduces on another machine at all is that almost everyone finally agreed on what a number is.

20.2 Machine Epsilon and the Precision Intrinsics

If the representable numbers form a grid, the natural next question is: how fine is the grid? How far apart are neighbouring doubles? The answer is the most important single quantity in numerical computing, and Fortran hands it to you as an intrinsic.

Definition (machine epsilon). Machine epsilon, written $\varepsilon_{\text{mach}}$, is the gap between $1.0$ and the next larger representable number. For double precision it is $2^{-52} \approx 2.22 \times 10^{-16}$; for single precision, $2^{-23} \approx 1.19 \times 10^{-7}$. Fortran returns it with the intrinsic epsilon(x), where the value of x is irrelevant — only its kind matters, so epsilon(1.0_dp) reports the double-precision epsilon. Intuitively, machine epsilon is the relative resolution of the format: it is the size of the smallest change you can make to a number near 1, and equivalently the reason a real(dp) carries about $-\log_{10}(\varepsilon_{\text{mach}}) \approx 16$ significant decimal digits.

A closely related quantity deserves its own name, because the two are often confused.

Definition (unit roundoff). The unit roundoff $u$ is the maximum relative error introduced by rounding a real number to the nearest representable value. Under the default round-to-nearest rule, a number is never more than half a gap from the nearest grid point, so $u = \varepsilon_{\text{mach}}/2 = 2^{-53} \approx 1.11 \times 10^{-16}$ for double precision. The rule to memorize: every basic floating-point operation returns the exactly-rounded result, so its relative error is at most $u$. That single bound, applied operation by operation, is how all round-off analysis begins — including the Project Checkpoint's error budget for the heat solver.

Fortran gives you a small family of numeric inquiry intrinsics that report the properties of a kind without your having to remember any of the numbers above. They are the professional's toolkit for reasoning about precision portably:

Here they are in action, with every printed value computed by hand:

program precision_intrinsics
  use kinds, only: dp
  implicit none

  print '(a, es23.15e3)', 'epsilon(1.0_dp)  = ', epsilon(1.0_dp)
  print '(a, es23.15e3)', 'huge(1.0_dp)     = ', huge(1.0_dp)
  print '(a, es23.15e3)', 'tiny(1.0_dp)     = ', tiny(1.0_dp)
  print '(a, es23.15e3)', 'spacing(1.0_dp)  = ', spacing(1.0_dp)
  print '(a, es23.15e3)', 'spacing(1.0e6_dp)= ', spacing(1.0e6_dp)
  print '(a, i0)',        'precision(1.0_dp)= ', precision(1.0_dp)
  print '(a, i0)',        'digits(1.0_dp)   = ', digits(1.0_dp)
end program precision_intrinsics

$ gfortran -std=f2018 -Wall precision_intrinsics.f90 -o pintr && ./pintr
epsilon(1.0_dp)  =  2.220446049250313E-016
huge(1.0_dp)     =  1.797693134862316E+308
tiny(1.0_dp)     =  2.225073858507201E-308
spacing(1.0_dp)  =  2.220446049250313E-016
spacing(1.0e6_dp)=  1.164153218269348E-010
precision(1.0_dp)= 15
digits(1.0_dp)   = 53

The two spacing lines carry the deepest idea in this section. spacing(1.0_dp) equals machine epsilon, $2^{-52}$ — the gap right above $1.0$. But spacing(1.0e6_dp) is $2^{-33} \approx 1.16 \times 10^{-10}$, about a million times larger. That is not a coincidence: it is exactly a million times larger, because $10^{6}$ is about a million times larger than $1$. The grid is not uniform. Near $1$, doubles are spaced about $2.2\times10^{-16}$ apart; near a million, about $1.2\times10^{-10}$ apart; near $10^{16}$, they are spaced more than $1$ apart, so consecutive integers can no longer all be represented. The spacing scales with the magnitude of the number — which is the whole point of floating point, and the source of its name: the decimal point floats to keep a constant number of significant digits, trading absolute precision for enormous range.

💡 Intuition — ULP, the unit in the last place. The spacing(x) intrinsic returns what numerical analysts call one ULP (unit in the last place): the value of a one-bit change in the significand at x's magnitude. Because the significand has a fixed 53 bits but the exponent slides, one ULP is a constant fraction of x (about $\varepsilon_{\text{mach}}$ of it) but a growing absolute amount as x grows. When we say a result is "accurate to a few ULPs," we mean it lands within a few grid points of the true answer — the natural unit of floating-point error, because it automatically scales with the numbers involved.

The nearest intrinsic lets you step along the grid one point at a time, and reveals a subtlety worth seeing. The gap above $1.0$ is $2^{-52}$, but the gap below $1.0$ is only $2^{-53}$ — half as large — because just below a power of two the exponent drops by one and the grid suddenly doubles in density:

program nearest_demo
  use kinds, only: dp
  implicit none
  print '(a, es13.6)', 'gap above 1.0 : ', nearest(1.0_dp,  1.0_dp) - 1.0_dp
  print '(a, es13.6)', 'gap below 1.0 : ', 1.0_dp - nearest(1.0_dp, -1.0_dp)
end program nearest_demo

$ gfortran -std=f2018 -Wall nearest_demo.f90 -o near && ./near
gap above 1.0 :  2.220446E-16
gap below 1.0 :  1.110223E-16

⚡ Performance Note — precision has a price, and it is paid in bandwidth. It is tempting to think "double is safer, so always use double." But a single-precision real is 4 bytes and a double is 8, and as Chapter 3 foreshadowed and Chapter 27 will measure, most serious numerical codes are memory-bandwidth-bound — they spend their time waiting for numbers to arrive from memory, not doing arithmetic. Halving the bytes can nearly double the speed of such a loop, and a SIMD vector register holds twice as many singles as doubles, doubling the arithmetic per instruction on top. Performance is not accidental: choosing precision is choosing speed, and §20.6 makes the trade explicit. The reason this book still defaults to dp is that scientific accuracy usually needs those extra nine decimal digits — but "usually" is a decision, not a reflex.

🔄 Check Your Understanding. 1. What does epsilon(1.0_dp) return, and why does the value of its argument not matter? 2. Is spacing(1.0e10_dp) larger or smaller than spacing(1.0_dp)? By roughly what factor? 3. Your code needs "about 12 significant digits." Is single precision enough? Is double?

Answers

(1) It returns $2^{-52} \approx 2.22\times10^{-16}$, the gap above $1.0$ in double precision. Only the kind of the argument matters — epsilon is a property of the type, not of any particular value — so epsilon(1.0_dp) and epsilon(3.7e9_dp) return the same number. (2) Larger, by a factor of about $10^{10}$: spacing scales with magnitude. (3) Single precision guarantees only ~6 digits (precision returns 6), nowhere near enough; double guarantees 15, comfortably enough. This is exactly the kind of decision §20.6 formalizes.

20.3 Catastrophic Cancellation, and How to Avoid It

Rounding a number to the grid costs you at most half a ULP — a relative error of $u \approx 10^{-16}$, which sounds harmless, and usually is. The danger is not a single rounding; it is an arrangement of arithmetic that amplifies those tiny relative errors into enormous ones. There is one operation that does this above all others, and it has a name.

Definition (catastrophic cancellation). Catastrophic cancellation is the severe loss of significant digits that occurs when you subtract two nearly equal floating-point numbers. The subtraction itself is often exact; the catastrophe is that it strips away the leading digits the two numbers agreed on, promoting their small, previously-negligible rounding errors into the leading digits of the result. If two numbers agree to 15 digits and you subtract them, the difference may have only 1 correct digit — you started with double precision and ended with barely one significant figure.

The clearest way to feel this is to work a small example in decimal with a fixed number of digits, the way a limited-precision machine does. Suppose we keep 8 significant decimal digits, and we want $x = 123456.79 - 123456.78 = 0.01$. Each operand is a perfectly good 8-digit number. But suppose each carries a rounding error in its last digit — each is really the true value $\pm$ a few units in that 8th place. Subtract them: the six leading digits 123456 cancel exactly, and what survives is $0.01 \pm (\text{those last-digit errors})$. The errors that were one part in $10^{8}$ of the operands are now a large fraction of the result. You went from 8 good digits to perhaps 1. No bits were lost in the subtraction; the information was never there — it was hiding behind digits that cancelled.

In real double precision the same thing happens far out in the sixteenth digit, invisibly, and two demonstrations make it concrete. Both compile and run:

program cancellation
  use kinds, only: dp
  implicit none
  real(dp), parameter :: big = 1.0e17_dp
  real(dp) :: tiny_bit, recovered

  ! (1) Absorption breaks associativity: 1.0 is lost beside 1.0e17.
  print '(a, f4.1)', '(big + 1) - big = ', (big + 1.0_dp) - big
  print '(a, f4.1)', 'big - big + 1   = ', (big - big) + 1.0_dp

  ! (2) A small quantity added to 1 and then subtracted away is gone.
  tiny_bit  = epsilon(1.0_dp) / 4.0_dp        ! 2^-54, well below half a ULP at 1.0
  recovered = (1.0_dp + tiny_bit) - 1.0_dp    ! should be tiny_bit; is exactly 0
  print '(a, es9.2)', 'tiny_bit        = ', tiny_bit
  print '(a, es9.2)', '(1+tiny) - 1    = ', recovered
end program cancellation

$ gfortran -std=f2018 -Wall cancellation.f90 -o cancel && ./cancel
(big + 1) - big =  0.0
big - big + 1   =  1.0
tiny_bit        =  5.55E-17
(1+tiny) - 1    =  0.00E+00

The first pair of lines proves floating-point addition is not associative: $(big + 1) - big$ gives 0.0 while $big - big + 1$ gives 1.0, from the same three numbers in a different order. The reason is absorption: $10^{17}$ has a spacing (ULP) larger than $1$, so big + 1.0 rounds right back to big, and subtracting big then yields zero. The 1 was destroyed the instant it was added to something so much larger. The second pair shows the same effect as pure cancellation: tiny_bit is a real, nonzero $2^{-54}$, but adding it to $1.0$ and subtracting $1.0$ back returns exactly 0 — the addition rounded it away, and the subtraction exposed the loss.

🐛 Find the Bug. This function is meant to compute $f(x) = \sqrt{x^2 + 1} - x$. It is correct for small x but returns nonsense — eventually a flat 0.0 — for large x. Why, and how would you fix it without changing what it computes?

fortran pure function f_naive(x) result(y) real(dp), intent(in) :: x real(dp) :: y y = sqrt(x*x + 1.0_dp) - x end function f_naive

Diagnosis and fix

For large x, $\sqrt{x^2 + 1}$ is only a whisper larger than x — they agree to nearly all 16 digits — so the subtraction is catastrophic cancellation, and for x beyond about $10^{8}$, $x^2 + 1$ rounds to $x^2$ (absorption), $\sqrt{x^2}$ is exactly x, and the result collapses to 0. The fix is algebra, not arithmetic: multiply by the conjugate, $\sqrt{x^2+1} - x = \dfrac{(\sqrt{x^2+1}-x)(\sqrt{x^2+1}+x)}{\sqrt{x^2+1}+x} = \dfrac{1}{\sqrt{x^2+1}+x}$. The rewritten form y = 1.0_dp / (sqrt(x*x + 1.0_dp) + x) computes the same mathematical value with an addition where the subtraction used to be — no cancellation, accurate for all x > 0. This conjugate trick is the single most useful cancellation cure you will learn.

The most famous casualty of cancellation is the quadratic formula. To solve $ax^2 + bx + c = 0$, the textbook writes $x = \dfrac{-b \pm \sqrt{b^2 - 4ac}}{2a}$. When $b^2 \gg 4ac$, the square root is very close to $|b|$, and one of the two roots is computed as a difference of nearly equal numbers — the small root is annihilated. The cure again is algebra: compute the large root first (the one where the $\pm$ becomes an addition, with no cancellation), then get the small root from the identity $x_1 x_2 = c/a$, so $x_{\text{small}} = c / (a\,x_{\text{large}})$. Same roots, no subtraction of near-equals.

💡 Intuition — you cannot recover information that cancelled. Catastrophic cancellation is not a rounding you can fix by rounding more carefully; the significant digits are simply gone, cancelled away with the digits the operands shared. The only real remedy is to never form the difference — to rearrange the mathematics so the near-equal subtraction disappears, as the conjugate trick and the quadratic-root swap both do. When you cannot avoid it, you must add precision (§20.6) so that the digits that survive are still enough. Recognizing cancellation before you code the formula, and reaching for the algebraic identity, is a defining skill of a computational scientist.

🔗 Connection. Cancellation is not an academic worry; it is why production numerical libraries look the way they do. LAPACK's routines (Chapter 21) are painstakingly arranged to avoid it, and the finite-difference formulas of Chapter 22 live in constant tension with it: a derivative approximation $\big(f(x+h) - f(x)\big)/h$ subtracts two nearly equal values, so shrinking h to reduce truncation error eventually increases rounding error as cancellation takes over — there is an optimal h, and it is set by machine epsilon.

20.4 Overflow, Underflow, NaN, and Inf

The grid of representable numbers is not only discrete; it is finite. There is a largest finite double, huge(1.0_dp) $\approx 1.8\times10^{308}$, and a smallest normal positive one, tiny(1.0_dp) $\approx 2.2\times10^{-308}$. Push a computation past those edges and IEEE 754 has defined, specific behaviours — and two special values that every scientific programmer must recognize on sight.

Definition (overflow and underflow). Overflow occurs when a result's magnitude exceeds the largest finite representable number; under IEEE 754 the result becomes a signed infinity. Underflow occurs when a nonzero result is smaller in magnitude than tiny — smaller than the smallest normal number; the result is represented, with reduced precision, as a subnormal number, and if it is smaller still it flushes to zero. Overflow loses the number entirely (it becomes Inf); underflow degrades gracefully toward zero.

Definition (subnormal). A subnormal (or denormal) number fills the gap between tiny and zero. Below $2^{-1022}$ the format gives up the hidden leading 1 and lets the significand's leading bits be zero, so it can represent still-smaller magnitudes — down to about $4.9\times10^{-324}$ — but with fewer and fewer significant bits. Subnormals keep underflow gradual instead of abrupt, at some cost in speed on hardware that handles them in microcode.

Definition (NaN and Inf). Inf (infinity) is the result of overflow or of a nonzero number divided by zero; it has a sign, obeys sensible rules ($\text{Inf} + 1 = \text{Inf}$, $1/\text{Inf} = 0$), and compares greater than every finite number. NaN ("Not a Number") is the result of a genuinely undefined operation — $0/0$, $\text{Inf} - \text{Inf}$, $\sqrt{-1}$ in real arithmetic, $\log(-1)$. Its defining property is contagion and incomparability: any arithmetic involving a NaN produces a NaN, and — the property that surprises everyone — a NaN is not equal to anything, including itself. NaN == NaN is .false.. That last fact is not a curiosity; it is the portable, standard way to detect a NaN.

Modern Fortran provides a standard, intrinsic module for working with these values safely, and its existence is a small proof that modern Fortran is a modern language — most languages leave you poking at bit patterns, while Fortran 2003 standardized the interface:

program special_values
  use kinds, only: dp
  use, intrinsic :: ieee_arithmetic, only: ieee_is_nan, ieee_is_finite
  implicit none
  real(dp) :: x, zero, inf, not_a_number, h

  x    = 1.0_dp
  zero = 0.0_dp
  inf          = x    / zero        ! 1.0 / 0.0  ->  +Inf
  not_a_number = zero / zero        ! 0.0 / 0.0  ->  NaN
  h            = huge(1.0_dp)       ! a variable, so h + h overflows at run time

  print '(a, l1)', 'ieee_is_finite(1.0_dp)       = ', ieee_is_finite(x)
  print '(a, l1)', 'ieee_is_finite(inf)          = ', ieee_is_finite(inf)
  print '(a, l1)', 'ieee_is_finite(h + h)        = ', ieee_is_finite(h + h)
  print '(a, l1)', 'ieee_is_nan(inf)             = ', ieee_is_nan(inf)
  print '(a, l1)', 'ieee_is_nan(0/0)             = ', ieee_is_nan(not_a_number)
  print '(a, l1)', 'NaN == NaN ?                 = ', (not_a_number == not_a_number)
  print '(a, l1)', 'inf > huge(1.0_dp) ?         = ', (inf > huge(1.0_dp))
end program special_values

$ gfortran -std=f2018 -Wall special_values.f90 -o special && ./special
ieee_is_finite(1.0_dp)       = T
ieee_is_finite(inf)          = F
ieee_is_finite(h + h)        = F
ieee_is_nan(inf)             = F
ieee_is_nan(0/0)             = T
NaN == NaN ?                 = F
inf > huge(1.0_dp) ?         = T

Study the fifth and sixth lines together. ieee_is_nan(not_a_number) is T — the reliable test — while not_a_number == not_a_number is F, the self-inequality that defines a NaN. Before ieee_arithmetic existed, programmers detected NaN precisely by writing if (x /= x), which is still a correct and portable idiom; the intrinsic simply says what you mean. Note also that h + h (twice the largest finite double) overflows to Inf, so ieee_is_finite reports F: overflow is not an error that stops your program by default — it silently poisons the result, which is exactly why you must go looking for it.

⚠️ Common Pitfall — NaN and Inf are silent by default. A division by zero or a $\sqrt{-1}$ does not crash your program or print a warning; it quietly produces Inf or NaN, which then spreads through every subsequent calculation (any arithmetic touching a NaN yields NaN) until, hours into a run, your output file is full of NaN and you have no idea where it started. The defensive habit: during development, ask the compiler to trap these events and stop the moment one occurs.

🔗 Connection. That trap is a compiler flag you already met in Chapter 13: -ffpe-trap=invalid,zero,overflow turns an invalid operation (which produces NaN), a division by zero (which produces Inf), and an overflow into an immediate crash with a backtrace pointing at the exact line, converting a silent poisoning into a loud, locatable failure. Build with it during development and remove it for production runs (where a stray underflow-to-zero may be harmless and you do not want the crash). Pairing -ffpe-trap with ieee_is_nan checks at the boundaries of your subroutines is the standard way scientific codes stay honest.

🔄 Check Your Understanding. 1. What does 0.0_dp / 0.0_dp produce? What about 1.0_dp / 0.0_dp? 2. Why is if (x == x) a valid — if cryptic — way to ask "is x not a NaN?" 3. Your long simulation ends with an output file full of NaN. What one compiler flag would have caught the first one at its source?

Answers

(1) 0.0/0.0 is NaN (an undefined form); 1.0/0.0 is +Inf (overflow-like, a defined signed infinity). (2) Because a NaN is the only value not equal to itself, x == x is true for every ordinary number and false only when x is a NaN — so it tests "not a NaN." Prefer .not. ieee_is_nan(x) for readability. (3) -ffpe-trap=invalid,zero,overflow (from Chapter 13), which halts with a backtrace at the first offending operation.

20.5 Numerical Stability and Conditioning

We now have the vocabulary for the deepest question in Part V: when can you trust a computed answer? The surprising truth is that the answer depends on two separate things — a property of the problem and a property of the algorithm — and confusing them is the most common mistake in numerical work.

Definition (conditioning). Conditioning measures how sensitive a problem's answer is to small changes in its inputs — it is a property of the mathematics, independent of any algorithm or computer. A problem is well-conditioned if small relative changes in the input cause only small relative changes in the output, and ill-conditioned if they cause large ones. The amplification factor is the condition number: a condition number of $10^{k}$ means you can expect to lose about $k$ significant digits just from the input's rounding error, no matter how perfectly you compute. Evaluating a function near a steep root, or solving a linear system whose rows are nearly parallel, are classically ill-conditioned.

Definition (numerical stability). Numerical stability measures how much error an algorithm adds on top of what the problem's conditioning already forces — it is a property of the method, not the problem. A stable algorithm produces an answer about as good as the conditioning allows; an unstable one manufactures extra error through its own arithmetic (typically by internal catastrophic cancellation or by amplifying round-off at each step), returning garbage even for a well-conditioned problem. The quadratic formula's naive form from §20.3 is unstable; the rearranged form is stable — for the very same, perfectly well-conditioned, problem.

The distinction is worth a sharp statement, because it governs how you diagnose a bad result:

🚪 Threshold Concept — conditioning is the problem's fault; stability is yours. If a computation gives a poor answer, there are two possible culprits and they call for opposite responses. If the problem is ill-conditioned (high condition number), no algorithm can save you — the accurate answer is genuinely not determined by the inputs you have, and you must reformulate the problem or gather better data. If the problem is well-conditioned but your algorithm is unstable, the fix is entirely in your hands: choose a better-arranged method, exactly as the conjugate trick and the quadratic-root swap did. The professional habit is to ask, in order: "Is my problem well-conditioned?" and only then "Is my algorithm stable?" A great deal of wasted effort — and many wrong papers — come from tuning a stable algorithm to fix what is really an ill-conditioned problem, or blaming the problem for what is really an unstable method.

This pairing is the through-line of the rest of Part V, and it is worth previewing where each idea returns so you can watch for it:

• Linear algebra (Chapter 21): the condition number of a matrix predicts how many digits you lose solving $A\mathbf{x} = \mathbf{b}$, and it is why you call LAPACK — its algorithms are provably stable — instead of coding textbook Gaussian elimination, which is not.
• Differentiation and integration (Chapter 22): the trade-off between truncation error (smaller step = better) and round-off from cancellation (smaller step = worse) sets an optimal step size, straight from this chapter's arithmetic.
• Differential equations (Chapter 23 and Chapter 24): a time-stepping scheme is stable only if each step does not amplify the previous step's error. For the heat solver this becomes a hard limit on the timestep — the CFL condition of Chapter 24 — and it is the same stability idea, now applied to marching forward in time rather than to a single formula.

💡 Intuition — error is a budget you spend once. Think of every computed result as starting with a fixed budget of about 16 significant digits (in double precision). Conditioning charges you $k$ digits up front, just for having the problem you have. Instability charges you more, digit by digit, as the algorithm runs. Accumulated round-off over many steps charges a little at each step. Your job as a computational scientist is to keep the total spend below your budget — by picking a well-conditioned formulation, a stable algorithm, and (the subject of §20.6) enough precision that the budget was large enough to begin with.

20.6 Choosing Precision Deliberately

We can now answer the question Chapter 3 posed and postponed: single, double, or quadruple — how do you choose? The answer is an engineering trade-off, and making it consciously is one of the things a computational scientist is paid to do.

The three axes are accuracy, memory, and speed, and they pull against each other:

Defining them portably follows the Chapter 3 pattern exactly — describe the precision you need and let the compiler supply the kind:

module precisions
  implicit none
  private
  public :: sp, dp, qp
  integer, parameter :: sp = selected_real_kind(6, 37)      ! single:  ~7 digits
  integer, parameter :: dp = selected_real_kind(15, 307)    ! double:  ~16 digits  (our workhorse)
  integer, parameter :: qp = selected_real_kind(33, 4931)   ! quad:    ~34 digits
end module precisions

Three practical warnings turn this table into judgement. First, quad is usually not hardware — on mainstream CPUs gfortran emulates real128 in software (via libquadmath), so it can be tens to hundreds of times slower than double. Reach for it to check a double result or to rescue a specific ill-conditioned step, not as a blanket safety margin. Second, precision is not accuracy: computing in double does not make single-precision input data any more accurate — it only stops your arithmetic from adding error, so matching your working precision to your data's real accuracy is the honest choice. Third, the danger of single precision is accumulation: one single-precision operation loses only ~7 digits' worth, but a million of them, each feeding the next, can erode the answer catastrophically — which is precisely the situation a time-stepping simulation is in, and precisely why the heat solver will be dp.

🐍 Python Comparison. Python and NumPy make the same choice you do, usually without telling you. A bare Python float is always IEEE double — there is no single-precision Python scalar — so pure-Python numerics are double whether you want it or not. NumPy, being built for exactly this trade, lets you pick: np.float32, np.float64 (the default), and np.float128 (which, like Fortran's quad, is often software-emulated and platform-dependent). When you wrap a Fortran kernel for Python with f2py (Chapter 15), the NumPy dtype on the Python side and the real kind on the Fortran side must agree — a real(dp) array pairs with float64, a real(sp) array with float32 — or the bytes are silently misread. Getting that correspondence right is the most common f2py mistake, and this chapter is why it matters.

🔄 Check Your Understanding. 1. Why is single precision often faster than double for a large array computation, despite doing the "same" arithmetic? 2. You need a reference answer accurate to 25 digits to check a double-precision code. Which kind, and what is the catch? 3. True or false: computing in double precision makes your single-precision measurements more accurate.

Answers

(1) Because single-precision reals are half the bytes, so a bandwidth-bound loop moves half as much data and a SIMD register holds twice as many values — often nearly a 2× speedup. (2) Quadruple precision (qp, ~34 digits); the catch is that it is usually software-emulated and can be 10–100× slower than double, so use it only for the check, not the production run. (3) False — precision limits the error your arithmetic adds; it cannot improve the accuracy of the input data itself.

Project Checkpoint

Every checkpoint so far has been about building the heat solver; this one is about justifying a decision inside it. Back in Chapter 3 you wrote integer, parameter :: dp = selected_real_kind(15, 307) into kinds.f90 and declared the plate's temperature field real(dp). At the time that was house style, adopted on faith. Now you can defend it — and defending your precision choice with an error budget is exactly the kind of reasoning Chapter 24's time-stepping will demand.

Here is the argument. The explicit heat solver marches the temperature field forward one small timestep at a time, and a realistic run takes many steps — call it $N = 10^{6}$. Each step performs a handful of floating-point operations per grid point, and each operation, by the unit-roundoff rule of §20.2, can add a relative error of up to $u = \varepsilon_{\text{mach}}/2$. In the worst case these round-off errors accumulate roughly linearly, so after $N$ steps the accumulated relative error is on the order of $N \cdot u$. That single expression decides the precision:

program checkpoint_precision
  use kinds, only: dp
  implicit none
  integer,  parameter :: n_steps  = 1000000                    ! a long run: 10^6 steps
  real(dp), parameter :: u_double = epsilon(1.0_dp) / 2.0_dp   ! unit roundoff, double
  real(dp), parameter :: u_single = 2.0_dp**(-24)              ! unit roundoff, single
  real(dp) :: err_double, err_single

  err_double = real(n_steps, dp) * u_double                    ! worst-case, double
  err_single = real(n_steps, dp) * u_single                    ! worst-case, single

  print '(a, i0)',     'time steps                 : ', n_steps
  print '(a, es10.3)', 'unit roundoff  (double)    : ', u_double
  print '(a, es10.3)', 'unit roundoff  (single)    : ', u_single
  print '(a, es10.3)', 'accumulated rel. err double: ', err_double
  print '(a, es10.3)', 'accumulated rel. err single: ', err_single
end program checkpoint_precision

$ gfortran -std=f2018 -Wall checkpoint_precision.f90 -o ckpt && ./ckpt
time steps                 : 1000000
unit roundoff  (double)    :  1.110E-16
unit roundoff  (single)    :  5.960E-08
accumulated rel. err double:  1.110E-10
accumulated rel. err single:  5.960E-02

The two bottom lines settle the question. In double precision, a million timesteps accumulate a worst-case relative error of about $1.1 \times 10^{-10}$ — ten digits still clean, far below any physical uncertainty in the problem; the round-off is invisible. In single precision the same run accumulates about $6 \times 10^{-2}$ — a 6% error, which would corrupt the second significant figure and render the simulation worthless. (These are pessimistic worst-case bounds; real errors, partly cancelling, typically grow like $\sqrt{N}\,u$ and are smaller — but the single-precision margin is so thin that even the optimistic estimate is uncomfortable.) The decision is not close: the solver must be dp. You have now bounded your solver's round-off and shown, with numbers, that the precision you chose on faith in Chapter 3 was the right one — and you have a template for making the same argument for any long accumulation you ever write. Keep kinds.f90 exactly as it is; it survives to the Chapter 38 capstone unchanged.

Summary

This chapter replaced the intuition that computer numbers behave like real numbers with an accurate model of what they actually are, and the tools to reason about the difference.

The three things to memorize. First, a real(dp) is a grid point, not a real number: epsilon(1.0_dp) $\approx 2.22\times10^{-16}$ is the grid's relative resolution, and you never test floats with ==. Second, subtracting nearly equal numbers is catastrophic — rearrange the algebra to avoid it. Third, NaN /= NaN; detect special values with ieee_arithmetic and hunt them with -ffpe-trap.

Spaced Review

Four questions revisiting Chapter 3 and Chapter 5. Answer from memory before opening the details.

1. (Ch. 3) Write the one line that portably defines the double-precision kind used throughout this book, and explain what each of its two arguments requests.
   

   Answer

   `integer, parameter :: dp = selected_real_kind(15, 307)`. The `15` requests at least 15 significant *decimal digits* of precision (satisfied by binary64's 53-bit significand); the `307` requests a decimal exponent range out to at least $10^{307}$ (satisfied by the 11-bit exponent). Describing the requirement, rather than hardcoding a kind number like `real(8)`, is what makes it portable.

2. (Ch. 3) Why does the literal 0.1 in real(dp) :: x = 0.1 still fail to store exactly one-tenth, even though you used double precision?
   

   Answer

   Precision (double) controls *how many* bits the significand has, not *which* numbers are representable. One-tenth is a repeating fraction in base 2, so it cannot be represented exactly in any finite binary significand — double simply rounds it to a closer grid point than single would. (A subtler point: writing `0.1` without a `_dp` suffix rounds it first in *single* precision, an even coarser approximation — always write `0.1_dp`.)

3. (Ch. 3) What is the difference between precision(x) and epsilon(x), and what does each return for real(dp)?
   

   Answer

   `precision(x)` returns the number of guaranteed decimal *digits* (`15` for `dp`); `epsilon(x)` returns machine epsilon, the *gap* above 1.0 as a real number ($2^{-52}\approx 2.22\times10^{-16}$ for `dp`). One counts digits, the other measures the smallest relative step — two views of the same 53-bit significand.

4. (Ch. 5) The heat solver stores temperature as a 2D real(dp), allocatable :: u(:,:) array and sums or sweeps over it. Why can the order in which you add its elements change the result, and which recurring theme does that illustrate?
   

   Answer

   Because floating-point addition is not associative (§20.3): summing a large array left-to-right accumulates differently than pairwise or in a different order, since each partial sum rounds. It illustrates the threshold concept that arrays hold grid points, not real numbers — and it is why reproducible parallel reductions ([Chapter 33](../../part-08-parallel-programming/chapter-33-openmp/index.md)) take care about summation order. Arrays are Fortran's superpower, but they obey floating-point arithmetic's laws, not grade-school algebra's.

What's Next

You now know what your numbers are, how much you can trust them, and how to keep an algorithm from squandering that trust. That foundation is exactly what the rest of Part V builds on. Next, Chapter 21 puts it to work on the first real numerical task: solving systems of linear equations. You will write $A\mathbf{x} = \mathbf{b}$ and discover why — armed with everything this chapter taught about conditioning and stability — you should almost never solve it with your own code, and should instead call LAPACK, the battle-tested Fortran library that sits underneath NumPy, MATLAB, and R. The condition number you met in §20.5 becomes a concrete diagnostic, and the precision you chose here becomes the precision LAPACK computes in. Let's do some linear algebra.