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.

Tags: Arbitrary-Precision Arithmetic, Data.Ratio, Haskell, Heron’s Method, Integer Arithmetic, Newton’s Method, Numerical Methods, Rational Approximations, Rational Arithmetic

This entry was posted on November 22, 2012 at 18:21 and is filed under Haskell, Numerical Methods. You can follow any responses to this entry through the RSS 2.0 feed. You can leave a response, or trackback from your own site.

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

Simon Read Says:
November 23, 2012 at 12:56 | Reply
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:
  November 25, 2012 at 11:04 | Reply
  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.
Simon Read Says:
November 24, 2012 at 08:59 | Reply
I have a few other mathematics ideas we might like to discuss, but where do we put them?

1) How do you know Newton-Raphson is converging? How do you MAKE it converge? I have discovered a test which makes it converge.

2) Power series for pi. I have discovered a very simple two-stage method which is better than any power series, (Ramanujan eat your heart out!) but not as fast as more modern quadratically-convergent methods. It might spark ideas on two-stage methods of computing functions.

3) Two-stage method of computing sine and cosine. The extra stage of computation gives you another degree of freedom which might let you optimise the accuracy/effort balance.

4) Puzzle: a goat in a circular field. Apparently this was part of a Cambridge University entrance exam for mathematics. It’s still a fun puzzle to solve.

Rod Carvalho

Heron’s method using integer arithmetic

Like this:

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

Leave a Reply Cancel reply

Follow “Rod Carvalho”