## Heron’s algorithm

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.

heron(x, ε, e) ⇒ e if abs(x - e2) < ε

heron(x, ε, e) ⇒ heron(x, ε, 1/2 × (e + (x ÷ e)))

heron(x:ℝnn, ε:ℝp) : ℝnn

heron(x, ε) ⇒ heron(x, ε, 1)

#### 1.1` `Tests

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 ✓

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:

Operators:

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

heron(x:FP64, ε:FP64) : 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)

#### 2.1` `Tests

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.