Heron's Method

The Idea

Heron of Alexandria described a method for computing $\sqrt{S}$ in the first century ^[1].

The iteration is self-correcting: overshoot and undershoot cancel each other, pulling every successive estimate toward $\sqrt{S}$ from above.

Why It Works — The Newton-Raphson Connection

Finding $\sqrt{S}$ is equivalent to finding the positive root of $f(x) = x^2 - S$. Newton-Raphson ^[2] linearises $f$ at the current guess and solves for the zero of the tangent line:

Substituting $f(x) = x^2 - S$ and $f'(x) = 2x$ gives exactly Heron's formula:

Convergence

Newton-Raphson converges quadratically when the initial guess is close enough to the root. Letting $\varepsilon_n = x_n - \sqrt{S}$ be the error at step $n$, a Taylor expansion of $f$ around $\sqrt{S}$ gives the error recurrence:

The number of correct decimal digits doubles with each iteration. The plot below shows this collapse — note how the curve drops off a cliff rather than declining steadily.

Worked Example — $\sqrt{2}$

Start with $x_0 = 1$. Applying the recurrence with $S = 2$:

After only three steps the result agrees with $\sqrt{2} = 1.41421356\ldots$ to eight significant figures. The 3D plot below shows why the iteration works geometrically — the surface $h(x, S) = (x + S/x)/2$ and the truth surface $z = \sqrt{S}$ intersect exactly along the fixed-point curve $x = \sqrt{S}$.

Implementation

The loop terminates when consecutive iterates differ by less than a chosen tolerance. A safe initial guess is $x_0 = S$ (large $S$) or $x_0 = 1$ (small $S$).

def heron_sqrt(s: float, tolerance: float = 1e-10) -> float:
    if s < 0:
        raise ValueError("Cannot take sqrt of a negative number")
    if s == 0:
        return 0.0
    x = s if s >= 1 else 1.0
    while True:
        x_next = (x + s / x) / 2
        if abs(x_next - x) < tolerance:
            return x_next
        x = x_next