Case Study 1: Reading and Porting a Convergence Loop

"An algorithm must be seen to be believed." — Donald Knuth

Executive Summary

Iterating until a quantity converges is the single most common loop in numerical computing — it is how you solve equations you cannot solve directly, and it is exactly the shape your heat solver's time loop will take. In this case study you are handed a small, correct piece of pure Python that computes a square root by the ancient Babylonian (Newton) method, and your job is to read it, understand its control flow, and port it faithfully to modern Fortran. Along the way you will confront the one decision that makes or breaks an iterative loop — how to write the exit condition — and see why the obvious choice (stop when nothing changes) is a floating-point trap that §4.1 warned you about.

Skills applied: the infinite do + exit loop (§4.3); relational and logical conditions (§4.1); the real-equality pitfall and tolerance-based tests (§4.1, forward to Chapter 20); translating control flow between Python and Fortran (§4.7); hand-tracing a loop to verify its output.

Background

The Babylonian method for $\sqrt{a}$ is beautifully simple: guess a value $x$, then repeatedly replace it with the average of $x$ and $a/x$. If the guess is too big, $a/x$ is too small, and their average is closer; if too small, the reverse. Each step roughly doubles the number of correct digits — this is Newton's method applied to $x^2 - a = 0$, and it converges quadratically. The update is:

$$ x_{\text{new}} = \tfrac{1}{2}\left(x + \frac{a}{x}\right) $$

You are given this working pure-Python implementation and asked to bring it into a Fortran code base:

def babylonian_sqrt(a, tol=1e-8):
    x = 1.0
    iters = 0
    while True:
        x_new = 0.5 * (x + a / x)
        iters += 1
        if abs(x_new - x) < tol:
            x = x_new
            break
        x = x_new
    return x, iters

It runs, and it is correct. The task is not to improve it but to reproduce it exactly in Fortran — same algorithm, same result — which forces you to read its control flow precisely.

Phase 1 — Read the Control Flow

Before writing a line of Fortran, name every control-flow element in the Python. There are only three, and they are all in this chapter:

The shape is the classic infinite loop with a mid-body exit: you cannot test the stopping condition at the top, because you do not have x_new until you have computed it, so the test sits in the middle of the body — compute, then check, then either exit or continue. This is precisely the pattern §4.3 introduced, and it is why Fortran's bare do … exit … end do exists.

Note the careful ordering: the code computes x_new, counts the iteration, tests, and only then overwrites x. The overwrite happens on both paths (inside the if before break, and after it), which keeps x holding the best estimate whether the loop exits or continues. A faithful port must preserve that ordering exactly.

Phase 2 — The Port

The translation is nearly mechanical once the control flow is named — the value of Phase 1. Here it is, computing $\sqrt{2}$, in the book's modern style: implicit none, real(dp), kind-suffixed literals.

program babylonian
  use, intrinsic :: iso_fortran_env, only: dp => real64
  implicit none
  real(dp) :: a, x, x_new, tol
  integer  :: iters

  a     = 2.0_dp
  tol   = 1.0e-8_dp
  x     = 1.0_dp
  iters = 0
  do
    x_new = 0.5_dp * (x + a / x)
    iters = iters + 1
    if (abs(x_new - x) < tol) then
      x = x_new
      exit
    end if
    x = x_new
  end do

  print '(a, i0)',    'iterations : ', iters
  print '(a, f12.8)', 'sqrt(2)    ~ ', x
end program babylonian

$ gfortran -std=f2018 -Wall -O2 babylonian.f90 -o babylonian && ./babylonian
iterations : 5
sqrt(2)    ~   1.41421356

The while True became a bare do; the break became exit; the condition carried over unchanged but for Fortran's .and.-family spelling (not needed here) and the mandatory then/end if. Every 1.0 became 1.0_dp so the arithmetic runs in the double precision the dp kind names — a habit from Chapter 3 that the Python, with its single float type, never had to think about.

Phase 3 — Trace the Convergence by Hand

The book's discipline is to never trust output we have not reasoned about, so trace the loop. Starting from $x = 1$, each pass computes $x_{\text{new}} = \tfrac{1}{2}(x + 2/x)$ and checks $|x_{\text{new}} - x|$ against $10^{-8}$:

The quadratic convergence is visible in the last column: the gap shrinks roughly as the square of the previous gap ($10^{-1} \to 10^{-2} \to 10^{-3} \to 10^{-6} \to 10^{-12}$), so the loop needs only 5 iterations. That matches the printed iterations : 5, and the final value rounds to 1.41421356 at eight decimals — the correctly-rounded double-precision $\sqrt{2}$.

Sanity check. A converged square root must satisfy $x^2 \approx a$. Squaring the result gives $1.41421356^2 \approx 2.0000000$, off from 2 by about $10^{-15}$ — the size of a rounding error in double precision, which is exactly what "converged" should mean.

Phase 4 — The Trap the Port Exposes

Here is where porting teaches you something. A tempting way to write the exit condition is "stop when the value stops changing":

if (x_new == x) exit      ! DANGEROUS as a general convergence test

For this problem it happens to work — Newton's method reaches an exact fixed point in double precision (pass 6 would produce a bit-identical x_new), so equality eventually holds. But as a general habit it is a bug waiting for a different function, where the iterate can oscillate in the last bit forever — flipping between two adjacent representable values, never exactly equal, and your loop never exits. This is the real-equality pitfall of §4.1 in its most dangerous form, because the loop looks correct and hangs only sometimes.

The robust condition is the one the Python used and the port kept: stop when the change falls below a tolerance, abs(x_new - x) < tol. The tolerance says "close enough," which is the only thing floating-point arithmetic can ever promise. Choosing tol well — not so tight the loop never converges, not so loose the answer is inaccurate — is a real skill, and its home is Chapter 20. For now, the lesson is narrow and permanent: an iterative loop's exit condition is a tolerance test, never an exact-equality test on reals.

Phase 5 — Why This Is the Heat Solver in Miniature

Step back and look at the shape you just ported: initialize a state, loop, update the state, test a stopping condition, exit. That is exactly the control flow of the time-stepping skeleton you build in this chapter's Project Checkpoint, and of the finished solver in Chapter 24 — the only differences are that the "state" is a whole temperature grid rather than a single number, and the stopping condition is "the plate has reached steady state" rather than "the root has converged." Master the loop here, on one number you can trace by hand, and the version with a million numbers holds no surprises.

Discussion Questions

1. The exit test sits in the middle of the loop body, not at the top. Why can this convergence loop not be written as a do while with the condition at the top? (Hint: what do you not yet have on the first pass?)
2. The Python uses while True: … break; the Fortran uses a bare do … exit. Are these the same construct with different spelling, or is there a semantic difference? Defend your answer.
3. Suppose you replaced tol = 1.0e-8_dp with tol = 1.0e-20_dp. Predict what happens to the iteration count, and whether the loop is guaranteed to terminate. (Consider the smallest gap double precision can represent.)

Your Turn: Extensions

• Option A. Add a safety cap: change the bare do to do while (iters < max_iters) or add if (iters > max_iters) exit inside, so a non-converging input cannot hang the program. What should the code print if it hits the cap? (This is the defensive pattern of Chapter 13.)
• Option B. Generalize the program to read a from the user and compute $\sqrt{a}$ for several values in a counted do loop, guarding against a < 0 with an if before you iterate.
• Option C. Port the cube root by the analogous Newton update $x_{\text{new}} = \tfrac{1}{3}!\left(2x + a/x^2\right)$, trace its convergence for $a = 27$ (the answer is 3), and confirm it still needs only a handful of iterations.

Key Takeaways

• The workhorse numerical loop is iterate until converged: an infinite do with a mid-body exit, because the stopping value does not exist until the body has run once.
• Porting forces you to name control flow precisely. while True → bare do; break → exit; the mid-body if carries over directly.
• The exit condition of an iterative loop is a tolerance test (abs(Δ) < tol), never an exact == on reals — the latter can hang forever on a value that oscillates in its last bit.
• This tiny loop is the heat solver's time loop in miniature: initialize, loop, update, test, exit. The structure is identical whether the state is one number or a million.