Exercises: Floating-Point Arithmetic

These exercises build the single most transferable numerical skill in the book: reasoning quantitatively about the error in a computation before it bites you. Some are pencil-and-paper; several ask you to predict a program's output and then compile it to check yourself — the fastest way to internalize how floating point really behaves. Reach for the precision intrinsics (epsilon, spacing, huge, tiny) rather than memorizing numbers.

Difficulty: ⭐ warm-up · ⭐⭐ standard · ⭐⭐⭐ deeper. Solutions: worked solutions to the daggered (†) and odd-numbered problems are in appendices/answers-to-selected.md; the computational ones are in code/exercise-solutions.f90. Try every problem before you look, and predict every output before you compile.

Part A — Warm-ups ⭐

20.1 † State how IEEE 754 double precision divides its 64 bits among sign, exponent, and significand, and say roughly how many significant decimal digits the significand buys you.

20.2 What is machine epsilon for double precision, and what does epsilon(1.0_dp) return? Why does the value of the argument to epsilon not matter?

20.3 † Explain, in terms of binary fractions, why 0.1 cannot be stored exactly as a real(dp). Name one other short decimal that also cannot, and one that can.

20.4 Distinguish overflow from underflow. What special value does each typically produce in IEEE arithmetic?

20.5 † Why is if (x == 0.3_dp) almost always a bug? Write the line you should use instead, and say what you would choose for the tolerance.

Part B — Type, Compile, and Run ⭐⭐

Predict the output of each, in writing, then compile and check. The gap between your prediction and the result is where the learning is.

20.6 † print '(l1)', (0.1_dp + 0.2_dp == 0.3_dp) — what prints, and why?

20.7 print '(f4.1)', (1.0e16_dp + 1.0_dp) - 1.0e16_dp — predict the value and explain it using the ULP at $10^{16}$.

20.8 † For x = 0.0_dp/0.0_dp and y = 1.0_dp/0.0_dp, predict ieee_is_nan(x), ieee_is_nan(y), and ieee_is_finite(y). Which of x, y is a NaN, and which an Inf?

20.9 Given that spacing(1.0_dp) is $2^{-52}$, predict spacing(2.0_dp) and spacing(1024.0_dp) as powers of two, and explain the pattern.

20.10 † Is (0.1_dp + 0.2_dp) - 0.3_dp positive, negative, or zero? Give its order of magnitude and express it as a power of two.

Part C — Find the Bug ⭐⭐

Each snippet is legal Fortran that misbehaves. Diagnose it and give the fix.

20.11 † This loop is meant to stop when x reaches 1.0, but it runs forever. Why, and how do you fix it?

real(dp) :: x = 0.0_dp
do while (x /= 1.0_dp)
   x = x + 0.1_dp
end do

20.12 A finite-difference derivative (f(x+h) - f(x)) / h is computed with h = 1.0e-20_dp and returns exactly 0.0 for a smooth f. Explain, and suggest a better h.

20.13 † A convergence test never fires:

if (abs(temp_new - temp_old) == 0.0_dp) exit

Diagnose it and rewrite it so the iteration actually stops when it has converged.

20.14 A program sums a large array two ways — a simple do loop and a loop over a permuted index — and the two totals differ in their last few digits. Is this a bug? Justify your answer with the right vocabulary.

Part D — Port It ⭐⭐

Translate from another language, and note where Fortran says the same thing more precisely.

20.15 † A NumPy user inspects their type with np.finfo(np.float64), reading off .eps, .max, and .tiny. Write the Fortran intrinsics that return the same three quantities for real(dp).

20.16 Port this naive Python quadratic-root function to Fortran, then rewrite it to be numerically stable for b*b >> 4*a*c:

def roots(a, b, c):
    d = (b*b - 4*a*c) ** 0.5
    return ((-b + d) / (2*a), (-b - d) / (2*a))

20.17 † In MATLAB a user relies on eps, realmax, and realmin. Give the Fortran intrinsic equivalents for real(dp), and note the one subtle difference in what MATLAB's eps means versus Fortran's epsilon.

Part E — Cancellation and Stability ⭐⭐⭐

20.18 † Rewrite $f(x) = \sqrt{x+1} - \sqrt{x}$ so it does not suffer catastrophic cancellation for large x. Show the algebra, and explain which subtraction you eliminated.

20.19 The expression $(1 - \cos x)/x^2$ loses accuracy for small x. Explain why, and give a rearrangement (a trig identity or a series) that is accurate near x = 0.

20.20 † For $x^2 - 10^{8}x + 1 = 0$ (so $a=1$, $b=-10^{8}$, $c=1$), the naive quadratic formula computes one root accurately and one disastrously. Identify which root is at risk, give the stable formula using $x_1 x_2 = c/a$, and state the two roots to a few significant figures.

20.21 Consider evaluating $f(x) = (e^{x} - 1)/x$ near x = 0. Is the problem ill-conditioned, or is the naive algorithm unstable? Justify, and name the intrinsic Fortran provides to sidestep the issue (hint: expm1 is a C function — does Fortran have a standard analogue, and what would you do without it?).

Part F — Back of the Envelope ⭐⭐

Order-of-magnitude reasoning. Show your work; the exponent matters, the mantissa does not.

20.22 † A simulation performs $10^{9}$ dependent floating-point operations, each with relative error at most the unit roundoff $u = 2^{-53}$. Bound the worst-case accumulated relative error (errors adding coherently) and the typical error (a random walk, $\sqrt{N}\,u$). Which is realistic, and why?

20.23 You need 12 correct significant digits after $10^{6}$ accumulating steps. Using $u \cdot N$ as a worst-case bound, decide whether single precision suffices, and whether double does.

20.24 † Exactly how many representable doubles lie in the interval $[1.0, 2.0)$? Explain your count from the bit layout.

20.25 A field is a $10{,}000 \times 10{,}000$ array. How much memory does it occupy in single versus double precision, and why might halving it more than halve your run time for a bandwidth-bound loop?

Part G — Design It (Extend the Solver) ⭐⭐

20.26 † Add a guard to kinds.f90 that makes the build fail loudly if the requested kind is unavailable — i.e., if selected_real_kind(15, 307) returns a negative value. Sketch the code and say where the check belongs.

20.27 Design a safe_divide(numer, denom) function that returns a supplied default (or calls error stop, per Chapter 13) when abs(denom) < tiny(1.0_dp), instead of silently producing Inf. When would the solver want this?

20.28 † The explicit heat solver's stable timestep scales as $dt \sim dx^2/(4\alpha)$ (you will derive the constant in Chapter 24). For $dx = 10^{-3}$ and $\alpha = 10^{-4}$, estimate dt, the number of steps in a 1-second simulation, and the worst-case accumulated round-off in double precision. Should you worry?

Part H — Interleaved ⭐⭐

Mixing this chapter with earlier ones.

20.29 † (Ch. 3) Explain why real(dp) :: x = 0.1 (note: no _dp suffix on the literal) is a double mistake — describe the two separate roundings that occur, and write the correct line.

20.30 (Ch. 5) The solver's field is real(dp), allocatable :: u(:,:). Explain why sum(u) may not equal a value you compute with your own nested do loop, and how the column-major layout of Chapter 5 relates to the order the elements are naturally visited.

20.31 † (Ch. 13) You suspect a 0.0/0.0 and a 1.0/0.0 are hiding somewhere in a long run. Which -ffpe-trap options would halt the program at the exact line of each, and why is that better than scanning the output for NaN afterward?

Solutions to the daggered and odd-numbered problems are in appendices/answers-to-selected.md; compilable solutions to the computational problems are in code/exercise-solutions.f90. Never trust a predicted output until you have compiled it — that is the whole discipline of this chapter.