Suppose we are given a positive rational number , whose square root we would like to compute. To be more precise, we would like to compute an approximation of .
Let be an approximation of . Hence, we have that is the geometric mean of and . If we replace the geometric mean of these two numbers by their arithmetic mean, we hopefully obtain a “better” approximation , as follows
If all goes well, then will be “sufficiently close” to for a “sufficiently large” integer . 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 and the sequence of approximations ?
Heron’s method using floating-point arithmetic
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
*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
Heron’s method using rational arithmetic
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
*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
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 can be expressed as a fraction , where and are integers and . Note that can (obviously?) be represented as a pair of integers . Since Haskell has arbitrary-precision integers, we can implement Heron’s method using pairs of arbitrary-precision integers to represent and the sequence of approximations .
Let and . The recurrence relation
can thus be rewritten in the following form
We then have two coupled recurrence relations
where , , and 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)
*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)
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 and by results in and being multiplied by . 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 , defined by . By construction, the zeros of function are and . Recall that Newton’s method uses the recurrence relation , where function is defined by
where is the first derivative of . Hence, we obtain the following
which is the recurrence relation used in Heron’s method.