## Heron’s method using integer arithmetic

Suppose we are given a positive rational number $y$, whose square root $\sqrt{y}$ we would like to compute. To be more precise, we would like to compute an approximation of $\sqrt{y}$.

Let $x_k \neq 0$ be an approximation of $\sqrt{y}$. Hence, we have that $\sqrt{y}$ is the geometric mean of $x_k$ and $y / x_k$. If we replace the geometric mean of these two numbers by their arithmetic mean, we hopefully obtain a “better” approximation $x_{k+1}$, as follows

$x_{k+1} = \displaystyle\frac{1}{2} \left( x_k + \frac{y}{x_k}\right)$

If all goes well, then $x_N$ will be “sufficiently close” to $\sqrt{y}$ for a “sufficiently large” integer $N$. This algorithm is the famous Heron’s method, devised by Heron of Alexandria some 2000 years ago. Some call it the Babylonian method, which suggests that the algorithm may be some 4000 years old. It’s not merely vintage, it’s antique.

Nonetheless, this is not a post on archaeology. We don’t want to sit around and talk about Heron’s method, we want to actually implement it. What programming language should we use? Let us use Haskell. How are we going to represent $y$ and the sequence of approximations $(x_k)_{k \in \mathbb{N}_0}$?

__________

Heron’s method using floating-point arithmetic

Though millions of people have already implemented Heron’s method using floating-point arithmetic (I implemented it in C in late 2000), one more implementation would do no harm:

y :: Floating a => a
y = 2

-- define Heron map
g :: Floating a => a -> a
g x = 0.5 * (x + y/x)

-- initial approximation
x0:: Floating a => a
x0 = 1

-- list of approximations
xs :: Floating a => [a]
xs = iterate g x0

Let us load this script into GHCi. Here’s a GHCi session:

*Main> take 6 xs
[1.0,1.5,1.4166666666666665,1.4142156862745097,
1.4142135623746899,1.414213562373095]
*Main> let x5 = xs !! 5
*Main> x5**2 - 2
-4.440892098500626e-16

We thus have the following approximation after five iterations

$x_5 = 1.414213562373095 \approx \sqrt{2}$

__________

Heron’s method using rational arithmetic

Let us now use rational numbers of arbitrary precision. The following script uses the Data.Ratio library:

import Data.Ratio

y :: Rational
y = 2

-- define Heron map
g :: Rational -> Rational
g x = 0.5 * (x + y/x)

-- initial approximation
x0:: Rational
x0 = 1 % 1

-- list of approximations
xs :: [Rational]
xs = iterate g x0

Let us load this script into GHCi. Here’s a GHCi session:

*Main> take 6 xs
[1 % 1,3 % 2,17 % 12,577 % 408,665857 % 470832,
886731088897 % 627013566048]
*Main> let x5 = xs !! 5
*Main> x5*x5 - 2
1 % 393146012008229658338304
*Main> 1 / 393146012008229658338304
2.5435842395854372e-24

Note that the error is eight orders of magnitude smaller than in the previous case (where we used finite-precision floating-point numbers). We pay for more accuracy with more computation. Thus, after five iterations, we have the following rational approximation

$x_5 = \displaystyle\frac{886731088897}{627013566048} \approx \sqrt{2}$

This is probably not news to you, my dear reader. Last April, I did, in fact, implement Heron’s method under the disguise of Newton’s method (see the appendix).

__________

Heron’s method using integer arithmetic

A rational number $y \in \mathbb{Q}$ can be expressed as a fraction $a / b$, where $a$ and $b$ are integers and $b \neq 0$. Note that $a / b$ can (obviously?) be represented as a pair of integers $(a,b)$. Since Haskell has arbitrary-precision integers, we can implement Heron’s method using pairs of arbitrary-precision integers to represent $y$ and the sequence of approximations $(x_k)_{k \in \mathbb{N}_0}$.

Let $y := a / b$ and $x_k := p_k / q_k$. The recurrence relation

$x_{k+1} = \displaystyle\frac{1}{2} \left( x_k + \frac{y}{x_k}\right)$

can thus be rewritten in the following form

$\displaystyle\frac{p_{k+1}}{q_{k+1}} = \frac{1}{2} \left( \frac{p_k}{q_k} + \frac{a}{b} \frac{q_k}{p_k} \right) = \frac{b \, p_k^2 + a \, q_k^2}{2 b \, p_k \, q_k}$

We then have two coupled recurrence relations

$\left[\begin{array}{c} p_{k+1}\\ q_{k+1}\end{array}\right] = \left[\begin{array}{c} b \, p_k^2 + a \, q_k^2\\ 2 b \, p_k \, q_k\end{array}\right]$

where $a$$b$$p_k$ and $q_k$ are integers. Here is a Haskell script that implements this vector recurrence relation:

a, b :: Integer
a = 2
b = 1

-- define Heron map
g :: (Integer,Integer) -> (Integer,Integer)
g (p,q) = (b * p^2 + a * q^2, 2 * b * p * q)

-- initial approximation
p0, q0 :: Integer
p0 = 1
q0 = 1

-- list of approximations
xs :: [(Integer,Integer)]
xs = iterate g (p0,q0)

Let us load this script into GHCi. Here’s a GHCi session:

*Main> take 6 xs
[(1,1),(3,2),(17,12),(577,408),(665857,470832),
(886731088897,627013566048)]
*Main> let x5 = xs !! 5
*Main> x5
(886731088897,627013566048)

Hence, after five iterations, we obtain the very same rational approximation we obtained in the previous case (where we used rational numbers of arbitrary-precision)

$x_5 = \displaystyle\frac{886731088897}{627013566048} \approx \sqrt{2}$

It should be noted, however, that this implementation using arbitrary-precision integers is different in one very important aspect from the one using arbitrary-precision rational numbers. What is the difference? Let GHCi answer this question:

*Main> import Data.Ratio
*Main Data.Ratio> 2 % 2
1 % 1
*Main Data.Ratio> 4 % 4
1 % 1
*Main Data.Ratio> 3 % 6
1 % 2

The main difference between the two implementations is that when we use arbitrary-precision rational numbers, there is automatic reduction of every fraction to its irreducible form, which requires computing the greatest common divisor (GCD) of the numerator and denominator of each rational number at each iteration. This has a computational cost, of course, but it also has the benefit of avoiding reducible fractions (which are inefficient representations of rational numbers). When we use pairs of arbitrary-precision integers, we do not compute the GCD (although we could, if we wanted to), but we may have to work with pairs of unnecessarily long integers (that correspond to reducible fractions), which is potentially dangerous.

For example, suppose that we replace the first three lines in the script above with the following three lines:

a, b :: Integer
a = 2000
b = 1000

We run the script again. Here’s yet another GHCi session:


*Main> take 6 xs
[(1,1),(3000,2000),(17000000000,12000000000),
(577000000000000000000000,408000000000000000000000),
(665857000000000000000000000000000000000000000000000,
470832000000000000000000000000000000000000000000000),
(8867310888970000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000000,
62701356604800000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000)]
*Main> let x5 = xs !! 5
*Main> x5
(8867310888970000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000000,
62701356604800000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000)

Frankly, I am shocked! Multiplying $a$ and $b$ by $1000$ results in $p_5$ and $q_5$ being multiplied by $10^{90}$. I may have made a mistake counting all those zeros, but it’s a lot more zeros than I expected.

__________

Appendix: from Heron to Newton

Interestingly, Heron’s method can be derived from the no less famous Newton’s method. Let us introduce a function $f : \mathbb{R} \to \mathbb{R}$, defined by $f (x) = x^2 - y$. By construction, the zeros of function $f$ are $\sqrt{y}$ and $-\sqrt{y}$. Recall that Newton’s method uses the recurrence relation $x_{k+1} = g ( x_k )$, where function $g : \mathbb{R} \to \mathbb{R}$ is defined by

$g (x) = x - \displaystyle\frac{f (x)}{f' (x)}$

where $f'$ is the first derivative of $f$. Hence, we obtain the following

$g (x) = x - \displaystyle\frac{f (x)}{f' (x)} = x - \displaystyle\left(\frac{x^2 - y}{2 x}\right) = \displaystyle\frac{1}{2} \left( x + \frac{y}{x}\right)$

which is the recurrence relation used in Heron’s method.

### 3 Responses to “Heron’s method using integer arithmetic”

In that case, you probably also know the very similar one for cube roots of $N$:

$x[i+1] = \displaystyle\frac{1}{3} \left( 2 x[i] + \frac{N}{ x[i]^2 } \right)$

which can also be derived from Newton-Raphson.

I’ve tried doing this in rational arithmetic and of course it works, but the numbers get rather out of hand: the number of digits multiplies by three each iteration, but the accuracy (number of accurate digits in floating point) only multiplies by two. That seems like the reward is not in proportion to the effort required.

• Rod Carvalho Says:

I’ve tried doing this in rational arithmetic and of course it works, but the numbers get rather out of hand (…) That seems like the reward is not in proportion to the effort required.

I have just written a post with some of my numerical experiments. I agree with you: the integers in the ratios just explode. It’s scary.

I suppose that the lesson to be learned from this is: with finite machines we can do finite-precision arithmetic comfortably, but doing arbitrary-precision arithmetic is a dangerous game.