# Analyzing Heron’s Algorithm

Finding the value of $\sqrt{a}$ for $a>0$ is equivalent to seeking a positive number $x$ whose square yields $a$: $x^2=a$. In other words, the solution to $x^2-a=0$.

We can find $\sqrt{4}$ by solving $x^2-4=0$. It is $2$ by mere inspection. $\sqrt{9}$ can also be obtained by solving $x^2-9=0$, again, by mere inspection. But ‘mere inspection’ can only go so far. For Example, what is $\sqrt{15241578750190521}$?

Method for finding the square root of a positive real number date back to the acient Greek and Babylonian eras. Heron’s algorithm, also known as the Babylonian method, is an algorithm named after the $1^{st}$ – century Greek mathematician Hero of Alexandria, It proceeds as follows:

1. Begin with an arbitrary positive starting value $x_0$.
2. Let $x_{k+1}$ be the average of $x_k$ and $\frac{a}{x_k}$
3. Repeat step 2 until the desired accuracy is achieved.

This post offers an analysis of Heron’s algorithm. Our aim is a better understanding of the algorithm through mathematics.

Let us begin by arbitrarily choosing a number $x_0>0$. If $x_0^2-a=0$, then $x_0$ is $\sqrt{a}$, and we have guessed the exact value of the square root. Otherwise, we are in one of the following cases:

Case 1: $x_0^2-a > 0 \implies x_0^2 > a \implies \sqrt{a} x_0 > a \implies \sqrt{a} > \frac{a}{x_0} \implies \frac{a}{x_0} < \sqrt{a} < x_0$

Case 2: $x_0^2-a <0 \implies x_0 < \sqrt{a} \implies \sqrt{a} x_0 < a \implies \sqrt{a} < \frac{a}{x_0} \implies x_0 < \sqrt{a} < \frac{a}{x_0}$

Both cases indicate that $\sqrt{a}$ lies somewhere between $\frac{a}{x_0}$ and $x_0$.

Let us define $e_0$ as the relative error of approximating $\sqrt{a}$ by $x_0$:

$e_0 =\frac{x_0-\sqrt{a}}{\sqrt{a}}$

The closer $e_0$ is to $0$, the better $x_0$ is as an approximation of $\sqrt{a}$.

Since $x_0>0, \sqrt{a}>0$,

$e_0 + 1 = \frac{x_0-\sqrt{a}}{\sqrt{a}} + 1 = \frac{x_0}{\sqrt{a}} > 0\quad\quad\quad(1)$

By (1),

$x_0 = \sqrt{a} (e_0+1)\quad\quad\quad(2)$

Let $x_1$ be the mid-point of $x_0$ and $\frac{a}{x_0}$:

$x_1 = \frac{1}{2} (x_0 + \frac{a}{x_0})\quad\quad\quad(3)$

and, $e_1$ the relative error of $x_1$ approximating $\sqrt{a}$,

$e_1 = \frac{x_1-\sqrt{a}}{\sqrt{a}}\quad\quad\quad(4)$

We have

$x_1-\sqrt{a}\overset{(3)}{=}\frac{1}{2}(x_0+\frac{a}{x_0}) - \sqrt{a} \overset{(2)}{=}\frac{1}{2}(\sqrt{a}(e_0+1)+\frac{a}{\sqrt{a}(e_0+1)})-\sqrt{a}=\frac{\sqrt{a}}{2}\frac{e_0^2}{e_0+1}\overset{(1)}{>}0\quad(5)$

$e_1 = \frac{x_1-\sqrt{a}}{\sqrt{a}} \overset{(5)}{=} \frac{\frac{\sqrt{a}}{2}\frac{e_0^2}{e_0+1}}{\sqrt{a}} = \frac{1}{2}\frac{e_0^2}{e_0+1}\overset{(1)}{ > }0\quad\quad\quad(6)$

By (5),

$x_1 > \sqrt{a} \implies \sqrt{a} x_1 > a \implies \sqrt{a} > \frac{a}{x_1}\implies \frac{\sqrt{a}}{x_1} < \sqrt{a} < x_1$

i.e., $\sqrt{a}$ lies between $\frac{a}{x_1}$ and $x_1$.

We can generate more values in stages by

$x_{k+1} = \frac{1}{2}(x_k + \frac{a}{x_k}), \;k=1, 2, 3, ...\quad\quad(8)$

Clearly,

$\forall k, x_k >0$.

Let

$e_k = \frac{x_k-\sqrt{a}}{\sqrt{a}}$.

We have

$x_k = \sqrt{a}(e_k+1)\quad\quad\quad(9)$

and

$e_{k+1} = \frac{x_{k+1}-\sqrt{a}}{\sqrt{a}}=\frac{\frac{1}{2}(\sqrt{a}(e_k+1) + \frac{a}{\sqrt{a}(e_k+1)})-\sqrt{a}}{\sqrt{a}}=\frac{1}{2}\frac{e_k^2}{e_k + 1}\quad\quad\quad(10)$

Consequently, we can prove that

$\forall k \ge 1, e_{k+1} > 0\quad\quad\quad(11)$

by induction:

When $k = 1, e_{1+1} \overset{(10)}{=} \frac{1}{2}\frac{e_1^2}{e_1+1} \overset{(6)} {>} 0$.

Assume when $k = p$,

$e_{p+1} =\frac{1}{2}\frac{e_p^2}{e_p^2+1} >0\quad\quad\quad(12)$

When $k = p + 1, e_{(p+1)+1} \overset{(10)}{=}\frac{1}{2}\frac{e_{p+1}^2}{e_{p+1} + 1}\overset{(12)}{>}0$.

Since (11) implies

$\forall k \ge 2, e_k > 0\quad\quad\quad(13)$

(10) and (13) together implies

$\forall k \ge 1, e_k > 0\quad\quad\quad(14)$

It follows that

$\forall k \ge 1, 0 \overset{(14)}{<} e_{k+1} \overset{(10)}{=} \frac{1}{2}\frac{e_k^2}{e_k+1} \overset{(14)}{<} \frac{1}{2}\frac{e_k^2}{e_k}=\frac{1}{2}e_k\quad\quad\quad(15)$

and

$\forall k \ge 1, 0 \overset{(14)}{<} e_{k+1} \overset{(10)}{=} \frac{1}{2}\frac{e_k^2}{e_k+1}\overset{(14)}{<}\frac{1}{2}\frac{e_k^2}{1}

(15) indicates that the error is cut at least in half at each stage after the first.

Therefore, we will reach a stage $n$ that

$0 < e_n < 1\quad\quad\quad(17)$

It can be shown that

$\forall k \ge 1, 0< e_{n+k} < e_n^{2^k}\quad\quad\quad(18)$

by induction:

When $k = 1$,

$0 \overset{(14)}{<} e_{n+1} \overset{(10)}{=}\frac{1}{2}\frac{e_n^2}{e_n+1} \overset{}{<} e_n^2=e_n^{2^1}\quad\quad\quad(19)$

Assume when k = p,

$0 < e_{n+p} < e_n^{2^p}\quad\quad\quad(20)$

When $k = p+1$,

$0 \overset{(14)}{<} e_{n+(p+1)} = e_{(n+p)+1} \overset{(16)}{<}e_{n+p}^2\overset{(20)}{<}(e_n^{2^p})^2 = e_n^{2^p\cdot 2} = e_n^{2^{p+1}}$

From (17) and (18), we see that eventually, the error decreases at an exponential rate.

Illustrated below is a script that implements the Heron’s Algorithm:

/*
a   - the positive number whose square root we are seeking
x0  - the initial guess
eps - the desired accuracy
*/
square_root(a, x0, eps) := block(
[xk],
numer:true,
xk:x0,
loop,
if abs(xk^2-a) < eps then return(xk),
xk:(xk+a/xk)/2,
go(loop)
)\$


Suppose we want to find $\sqrt{2}$ and start with $x_0 = 99$. Fig. 2 shows the Heron’s algorithm in action.

Fig. 1