Heron's algorithm for computing square roots

Heron’s algorithm, described in the 1st century CE by Heron of Alexandria, computes the square root of an input number iteratively, starting from an initial estimate, until the result is correct within a given tolerance. It is a special case of Newton’s method for finding roots of algebraic equations.

Uses Real numbers, Absolute value

Let heron(x:ℝ.nn, ε:ℝ.p, e:ℝ.nn) : ℝ.nn be the result of Heron’s algorithm for computing the square root of x up to tolerance ε, starting from estimate e.

The first step of the algorithm is to check if the current estimate is good enough, in which case it is the final result:
 heron[1]: heron(x:ℝ.nn, ε:ℝ.p, e:ℝ.nn) ⇒ e | abs(x - e2) < ε

Note that the tolerance applies to x and not to √(x).

Otherwise, a new estimate is computed by taking the average of e and x ÷ e:
 heron[2]: heron(x:ℝ.nn, ε:ℝ.p, e:ℝ.nn) ⇒ heron(x, ε, 1/2 × (e + (x ÷ e)))

For convenience, we also allow no initial estimate to be supplied, using a default value of 1:
 heron(x:ℝ.nn, ε:ℝ.p) : ℝ.nn
 heron(x:ℝ.nn, ε:ℝ.p) ⇒ heron(x, ε, 1)
The algorithm will always converge starting from 1, because √(x) is always in the interval [1..x] (for x > 1) or [x..1] (for x < 1). However, a better initial estimate makes it converge faster.

Example

The algorithm is straightforward to use:
 heron(2, 1/10) ⇒ 17/12

To see its convergence properties, it is more useful to look at the function
 convergence heronError(x:ℝ.nn) : ℝ.p → ℝ.nn,
(uses convergence Functions), computed as:
 convergence heronError(x:ℝ.nn)[ε:ℝ.p] ⇒ abs(heron(x, ε)2 - x)

A convergence plot for the square root of 2 (uses convergence Plotting):
 convergence linePlot(ε-values, heronError(2)[ε-values])
with
 convergence ε-values : Array1D(ℝ.p) and
 convergence ε-values ⇒ {[1, 1/2, 1/3, 1/5, 1/10, 1/100, 1/1000]}