Chapter 20 — Key Takeaways (Floating-Point Arithmetic)

A one-page reference for precision, rounding, and error. Pin it beside your solver.

The one idea

A real(dp) is not a real number — it is one point on a finite, unevenly spaced grid of about $2^{64}$ values. Every literal and every result is snapped to the nearest grid point. Arithmetic on the grid is commutative but not associative, most decimals are not representable, and == is a trap.

IEEE 754 formats

Format Fortran kind Bits (s / exp / signif.) Decimal digits epsilon Approx. max
single real32, selected_real_kind(6,37) 1 / 8 / 24 ~7 $2^{-23}\approx1.2\text{e-}7$ $3.4\text{e}38$
double (dp) real64, selected_real_kind(15,307) 1 / 11 / 53 ~16 $2^{-52}\approx2.2\text{e-}16$ $1.8\text{e}308$
quad real128, selected_real_kind(33,4931) 1 / 15 / 113 ~34 $2^{-112}\approx1.9\text{e-}34$ $1.2\text{e}4932$

Numbers worth memorizing (double precision)

  • epsilon(1.0_dp) $= 2^{-52} \approx 2.220446049250313\text{e-}16$ — the gap above 1.0.
  • unit roundoff $u = \varepsilon_{\text{mach}}/2 = 2^{-53} \approx 1.11\text{e-}16$ — max relative error of one operation.
  • huge(1.0_dp) $\approx 1.7977\text{e+}308$; tiny(1.0_dp) $\approx 2.2251\text{e-}308$.
  • precision(1.0_dp) = 15; digits(1.0_dp) = 53.
  • 0.1_dp + 0.2_dp $= 0.30000000000000004\ldots$ (not 0.3); the difference is $2^{-54}\approx5.55\text{e-}17$.

Precision intrinsics (reason about a kind portably)

Intrinsic Returns
epsilon(x) machine epsilon (gap above 1.0)
huge(x) / tiny(x) largest finite / smallest normal positive
spacing(x) one ULP at x (grows with magnitude)
nearest(x, s) next representable number toward sign(s)
precision(x) / digits(x) guaranteed decimal digits / significand bits

Rules that keep you out of trouble

  • Never compare floats with ==. Use abs(a - b) < tol, with tol scaled to the magnitudes.
  • Never write real literals without a kind suffix in numerical code: 0.1_dp, not 0.1.
  • Never subtract nearly equal numbers if you can rearrange the algebra to avoid it.
  • Absorption: a value below $\tfrac12\,\text{spacing}(S)$ added to a running sum $S$ is lost — so summation order matters.

Catastrophic cancellation — the cures

Unstable Stable rewrite
$\sqrt{x^2+1}-x$ $1/(\sqrt{x^2+1}+x)$ (conjugate)
$\sqrt{x+1}-\sqrt{x}$ $1/(\sqrt{x+1}+\sqrt{x})$
quadratic small root $\dfrac{-b+\sqrt{b^2-4ac}}{2a}$ large root first, then $x_{\text{small}} = c/(a\,x_{\text{large}})$
one-pass variance $\overline{x^2}-\bar x^2$ two-pass: subtract mean, then square

Special values

Situation Result Detect with
overflow, $1/0$ Inf (signed) .not. ieee_is_finite(x)
$0/0$, $\sqrt{-1}$, $\infty-\infty$ NaN ieee_is_nan(x) or x /= x

Key facts: NaN is contagious (any op with it → NaN) and NaN /= NaN. Inf > huge(1.0_dp) is .true.

Conditioning vs stability (§20.5)

  • Conditioning = the problem's sensitivity to input error (condition number $10^{k}$ → lose ~$k$ digits). You cannot beat it; reformulate or get better data.
  • Stability = extra error the algorithm adds. You can fix it — pick a better-arranged method.
  • Diagnose in order: "Is my problem well-conditioned?" then "Is my algorithm stable?"

Choosing precision (§20.6)

Want… Use Cost
speed/memory, ~6-digit data sp 1×
the scientific default dp 2×
a 34-digit reference oracle qp ~10–100× (usually emulated)

Danger of single: round-off accumulates over many dependent steps. Precision limits arithmetic error; it does not improve inaccurate input data.

Compile flags introduced

  • -ffpe-trap=invalid,zero,overflow (from Ch. 13) — crash with a backtrace at the first NaN/Inf.
  • esw.dEe edit descriptor (e.g. es23.15e3) — force a 3-digit exponent so E+308 prints with its E.

Project piece added this chapter

Not new code — a justified decision: an error budget ($N \cdot u$ over $\sim10^{6}$ steps) shows the heat solver's round-off is $\sim10^{-10}$ in dp (negligible) but $\sim6\%$ in single (ruinous). The kinds.f90 dp from Chapter 3 is confirmed as correct and stays unchanged to the Chapter 38 capstone.

Python bridge

Python float and NumPy float64 are IEEE binary64 — identical to real(dp). 0.1 + 0.2 misbehaves identically. When wrapping Fortran with f2py, real(dp) ↔ float64 and real(sp) ↔ float32 must match, or the bytes are misread.