Heron's algorithm

Heron’s algorithm

Konrad Hinsen

Heron’s algorithm computes the square root of an input number x iteratively, starting from an initial estimate e, until the result is correct within a given tolerance ε. It is a special case of Newton’s method for finding roots of algebraic equations.

Context heron
uses builtins/real-numbers

1 Heron’s algorithm using exact arithmetic

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 approximation is good enough, in which case it is the final result:

heron(x, ε, e) ⇒ e if 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(x, ε, e) ⇒ 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, ε) ⇒ heron(x, ε, 1)

The iteration starting from 1 will always converge but could well be inefficient.

1.1 Tests

We can use this algorithm with rational number arguments:

heron(2, 1/2)2 - 2 ⇒ 1/4 ✓
heron(2, 1/10)2 - 2 ⇒ 1/144 ✓
heron(2, 1/100)2 - 2 ⇒ 1/144 ✓
heron(2, 1/200)2 - 2 ⇒ 1/166464 ✓

We see that decreasing ε leads to better approximations of √2, which are always within the prescribed tolerance.

Context fp-heron
includes transformed heron

2 Heron’s algorithm using floating-point arithmetic

A floating-point version of Heron’s algorithm can be obtained by automatic conversion:

Sorts: FP64
Operators:

heron(x:FP64, ε:FP64, e:FP64) : FP64
heron(x:FP64, ε:FP64) : FP64

Vars: ε:FP64, x:FP64, e:FP64

Rules:

heron(x, ε, e) ⇒ e
if abs(x - e2) < ε
heron(x, ε, e) ⇒ heron(x, ε, 0.5 × (e + (x ÷ e)))
heron(x, ε) ⇒ heron(x, ε, 1.0)

Assets:

2.1 Tests

We can use this version with floating-point arguments:

heron(2.0, 0.5)2 - 2.0 ⇒ 0.25 ✓
heron(2.0, 0.5) - √(2.0) ⇒ 0.08578643762690485
abs(heron(2.0, 0.1)2 - 2.0) < 0.1 ⇒ true ✓
heron(2.0, 0.1) - √(2.0) ⇒ 0.002453104293571373
abs(heron(2.0, 0.01)2 - 2.0) < 0.01 ⇒ true ✓
heron(2.0, 0.01) - √(2.0) ⇒ 0.002453104293571373
abs(heron(2.0, 0.001)2 - 2.0) < 0.001 ⇒ true ✓
heron(2.0, 0.001) - √(2.0) ⇒ 2.123901414519125e-06

Again we see that decreasing ε leads to better approximations of √(2.0). The deviation is always smaller than the prescribed tolerance.

1	Heron’s algorithm using exact arithmetic
2	Heron’s algorithm using floating-point arithmetic