The Babylonian square root and bisection

Old Babylonian tablets give \(\sqrt{2}\) to about six decimal places. The rule Heron of Alexandria later wrote down averages a guess with \(a\) divided by the guess, \(x_{n+1} = \tfrac12(x_n + a/x_n)\). Bisection needs even less: it halves a bracket \([a, b]\) on which \(f\) changes sign, relying only on the intermediate value theorem. This script computes \(\sqrt{2}\) both ways and shows the bracket shrinking by exactly half per step, against the Babylonian rule’s doubling of correct digits.

import matplotlib.pyplot as plt
import numpy as np

from mathematicskit.numerical_analysis import Bisection, FixedPointIteration, estimate_convergence_order

root_exact = np.sqrt(2.0)

The Babylonian (Heron) iteration

\(g(x) = \tfrac12(x + 2/x)\) has \(\sqrt{2}\) as its fixed point. Starting from the tablet-friendly guess 1.5, six-place accuracy takes only three averaging steps.

babylonian = FixedPointIteration(lambda x: 0.5 * (x + 2.0 / x), x0=1.5, tol=1e-15).solve()
for n, x in enumerate(babylonian.history):
    print(f"Babylonian step {n}: x = {x:.15f}   error = {abs(x - root_exact):.1e}")

Babylonian step 0: x = 1.500000000000000   error = 8.6e-02
Babylonian step 1: x = 1.416666666666667   error = 2.5e-03
Babylonian step 2: x = 1.414215686274510   error = 2.1e-06
Babylonian step 3: x = 1.414213562374690   error = 1.6e-12
Babylonian step 4: x = 1.414213562373095   error = 2.2e-16
Babylonian step 5: x = 1.414213562373095   error = 2.2e-16

Bisection on \(f(x) = x^2 - 2\) over \([1, 2]\)

Each step keeps the half of the bracket where \(f\) changes sign, so the bracket width after \(n\) steps is exactly \(2^{-n}\).

bisection = Bisection(lambda x: x**2 - 2.0, 1.0, 2.0, tol=1e-13).solve()
print(f"bisection: root = {bisection.root:.13f} after {bisection.iterations} halvings")
print(f"Babylonian order ~ {estimate_convergence_order(babylonian.history[:-1], root_exact):.2f}")

bisection: root = 1.4142135623730 after 44 halvings
Babylonian order ~ 1.98

Bracket width vs. error

fig, ax = plt.subplots(figsize=(7, 4.5))
steps = np.arange(len(bisection.history))
bis_err = np.abs(bisection.history - root_exact)
bab_err = np.abs(babylonian.history - root_exact)
ax.semilogy(steps, 2.0 ** (-steps.astype(float)), "--", color="gray", label=r"bracket width $2^{-n}$")
ax.semilogy(steps, np.where(bis_err == 0, np.nan, bis_err), "o-", ms=3, color="darkorange", label="bisection midpoint")
ax.semilogy(np.arange(len(bab_err)), np.where(bab_err == 0, np.nan, bab_err), "s-", color="seagreen", label="Babylonian average")
ax.set_xlabel("step $n$")
ax.set_ylabel(r"$|x_n - \sqrt{2}|$")
ax.set_title(r"Computing $\sqrt{2}$: halving a bracket vs. Babylonian averaging")
ax.legend()
fig.tight_layout()

plt.show()