## Generalizing Heron’s method

Suppose we are given a positive rational number $y$. We would like to compute an approximation of the $n$-th root of $y$, which we denote by $\sqrt[n]{y}$. If $n = 2$, we can use Heron’s method to approximate $\sqrt{y}$, as we saw a few days ago. But, what if we have that $n > 2$? Is it possible to generalize Heron’s method?

Indeed, it is. Let us introduce a function $f_{n,y} : \mathbb{R} \to \mathbb{R}$, defined by $f_{n,y} (x) = x^n - y$, where $n \geq 2$ and $y > 0$. By construction, the positive zero of function $f_{n,y}$ is $\sqrt[n]{y}$. Do note that if $n$ is odd, then $-\sqrt[n]{y}$ is not a zero of $f_{n,y}$. Recall that Newton’s method uses the recurrence relation $x_{k+1} = g_{n,y} ( x_k )$, where function $g_{n,y} : \mathbb{R} \to \mathbb{R}$ is defined by

$g_{n,y} (x) = x - \displaystyle\frac{f_{n,y} (x)}{f_{n,y}' (x)}$

where $f_{n,y}'$ is the first derivative of $f_{n,y}$. Hence, we obtain

$g_{n,y} (x) = x - \displaystyle\frac{f_{n,y} (x)}{f_{n,y}' (x)} = x - \displaystyle\left(\frac{x^n - y}{n x^{n-1}}\right) = \displaystyle\frac{1}{n} \left( (n-1) x + \frac{y}{x^{n-1}}\right)$

Thus, we have the generalized Heron’s method

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

If we make $n = 2$, we obtain the standard Heron’s method, as we expected. The following Haskell script implements this generalized Heron’s method using arbitrary-precision rational numbers (included in the Data.Ratio library) to represent $y$ and the sequence of approximations $(x_k)_{k \in \mathbb{N}_0}$:

import Data.Ratio

-- define generalized Heron map
g :: Integer -> Rational -> Rational -> Rational
g n y x = (1 % n) * ((fromIntegral (n-1) * x) + (y / (x^(n-1))))

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

-- list of approximations of the n-th root of y
roots :: Integer -> Rational -> [Rational]
roots n y | n  > 1 = iterate (g n y) x0
| n <= 1 = error "Invalid value of n!!"

We load this script into GHCi. Let us now perform some numerical experiments.

__________

Experiment #1

We would like to compute a rational approximation of $\sqrt{2}$. We have $n = 2$ and $y = 2$. Here is a very brief GHCi session:

*Main Data.Ratio> take 6 (roots 2 2)
[1 % 1,3 % 2,17 % 12,577 % 408,665857 % 470832,
886731088897 % 627013566048]

We obtain the following rational approximation after five iterations

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

which is the same approximation I obtained twice before.

__________

Experiment #2

We would like to compute a rational approximation of $\sqrt{5}$. We have $n = 2$ and $y = 5$. Here is a very brief GHCi session:

*Main Data.Ratio> take 7 (roots 2 5)
[1 % 1,3 % 1,7 % 3,47 % 21,2207 % 987,4870847 % 2178309,
23725150497407 % 10610209857723]

We obtain the following rational approximation after six iterations

$x_6 = \displaystyle\frac{23725150497407}{10610209857723} \approx \sqrt{5}$

How good is this approximation? Let us check:

*Main Data.Ratio> (23725150497407 % 10610209857723)^2 - 5
4 % 112576553224922323902744729
*Main Data.Ratio> 4 / 112576553224922323902744729
3.5531377408652764e-26

Not bad! We could now obtain a rational approximation of the famous golden ratio $\varphi := \frac{1 + \sqrt{5}}{2}$, as follows:

*Main Data.Ratio> let xs = roots 2 5
*Main Data.Ratio> let phis = map (/2) (map (+1) xs)
*Main Data.Ratio> take 7 phis
[1 % 1,2 % 1,5 % 3,34 % 21,1597 % 987,3524578 % 2178309,
17167680177565 % 10610209857723]

We obtain the following rational approximation after six iterations

$\tilde{\varphi}_6 = \displaystyle\frac{17167680177565}{10610209857723} \approx \varphi$

How good is this approximation? Let us recall that the golden ratio is one of the solutions of the equation $x^2 - x - 1 = 0$. Hence, we can check how close to zero $\tilde{\varphi}_6^2 - \tilde{\varphi}_6 - 1$ is:

*Main Data.Ratio> let phi6 = 17167680177565 % 10610209857723
*Main Data.Ratio> phi6^2 - phi6 - 1
1 % 112576553224922323902744729
*Main Data.Ratio> 1 / 112576553224922323902744729
8.882844352163191e-27

I would say that is fairly close.

__________

Experiment #3

We now would like to compute a rational approximation of $\sqrt[3]{10}$. We have $n = 3$ and $y = 10$. Here is a very brief GHCi session:

*Main Data.Ratio> take 6 (roots 3 10)
[1 % 1,4 % 1,23 % 8,4909 % 2116,55223315303 % 25495981298,
83759169926117983945469262167029 % 38876457805393768546966848104041]

We obtain the following rational approximation after five iterations

$x_5 = \displaystyle\frac{83759169926117983945469262167029}{38876457805393768546966848104041} \approx \sqrt[3]{10}$

Is this approximation any good? Let us check:

*Main Data.Ratio> let x5 = (roots 3 10) !! 5
*Main Data.Ratio> x5^3
5876207108075025254603649658428809744645894398375582704
87986732898174481543266076676033790365389 % 58757060813
2677736272189042625165443146909261971794220978628909441
43919908404893015241356540921
*Main Data.Ratio> 5876207108075025254603649658428809744
6458943983755827048798673289817448154326607667603379036
5389 / 587570608132677736272189042625165443146909261971
79422097862890944143919908404893015241356540921
10.000852709004354

This is somewhat disappointing. The rational approximation $x_5$ is a ratio of two enormously long integers, but we do not have such a small error. Note that $x_5^3 - 10 \approx 10^{-3}$. We can do some quick error analysis. Let $f_{3,10} (x) = x^3 - 10$. The 1st order Taylor approximation of $f_{3,10}$ around $x_5$ evaluated at $x = \sqrt[3]{10}$ is

$f_{3,10} (\sqrt[3]{10}) \approx f_{3,10} (x_5) + f_{3,10}' (x_5) (\sqrt[3]{10} - x_5)$

and, since $f_{3,10} (\sqrt[3]{10}) = 0$, we obtain the estimate of the error

$x_5 - \sqrt[3]{10} \approx \displaystyle\frac{f_{3,10} (x_5)}{f_{3,10}' (x_5)}$

In other words, the output deviation $f_{3,10} (x_5) - f_{3,10} (\sqrt[3]{10})$ should be divided by the slope $f_{3,10}' (x_5)$ to obtain an estimate of the magnitude of the approximation error $x_5 - \sqrt[3]{10}$. Let us use GHCi yet once again:

*Main Data.Ratio> let error = (x5^3 - 10) / (3 * x5^2)
*Main Data.Ratio> error
1670089160826306272530773923851043922672595525468316978
5941152245094153072382174540074985393 % 272741620880842
8590518540423731512512814773109630109912907276686930689
12248623218095571481624481
*Main Data.Ratio> 1670089160826306272530773923851043922
6725955254683169785941152245094153072382174540074985393
/ 2727416208808428590518540423731512512814773109630109
91290727668693068912248623218095571481624481
6.123338108179484e-5

This is still quite a huge error, I would say.

__________

Bonus Experiment

Let us go crazy and attempt to compute a rational approximation of $\sqrt[100]{100}$. We have $n = 100$ and $y = 100$. We expect the rational numbers to quickly become ratios of astronomically long integers. Therefore, we will not be showing any rational approximations on GHCi. Here is something scary:

*Main Data.Ratio> let x3 = (roots 100 100) !! 3
*Main Data.Ratio> let error = (x3^100 - 100) / (100 * x3^99)
*Main Data.Ratio> let error_float = fromIntegral (numerator error)
/ fromIntegral (denominator error)
*Main Data.Ratio> error_float

After some 20 minutes without an output, I decided to abort this experiment. If you want to show x3, you can, but you will see thousands of digits. Thousands! You have been warned.

### 4 Responses to “Generalizing Heron’s method”

“rational approximation of \sqrt[3]{10}.
a very brief GHCi session:

*Main Data.Ratio> take 6 (roots 3 10)
1 % 1…

Note that x_5^3 – 10 \approx 10^{-3}.”

I see the problem: you started with 1/1 so used the first three iterations getting roughly near the root. If you start with (for example) 5/2, you’ll be much closer so allowing Newton-Raphson to do its best work. You’ll still have a lot of digits, though!

However, you can then use Haskell to get a rational approximation to your (super-multi-digit) approximation, so you might get something of similar accuracy but with fewer digits.

2. John Baez Says:

Are 17167680177565 and 10610209857723 Fibonacci numbers? The best rational approximations to the golden ratio are ratios of Fibonacci numbers.

• Rod Carvalho Says:

Prof. Baez,

Thanks for the insightful comment. I had not thought of using the Fibonacci sequence to compute the golden ratio, though now that you mention it, it seems obvious.

Indeed, those are Fibonacci numbers! Here’s the proof:

Prelude> let fibs = 1 : 1 : zipWith (+) fibs (tail fibs)
Prelude> take 20 fibs
[1,1,2,3,5,8,13,21,34,55,89,144,233,
377,610,987,1597,2584,4181,6765]
Prelude> 17167680177565 elem fibs
True
Prelude> 10610209857723 elem fibs
True

Surely they are fibonacci numbers. If you use rational arithmetic with Newton-Raphson to solve the polynomial $x^2-x-1=0$, then you get a selected tour of fractions using fibonacci numbers, *PROVIDED* you start with one pair in the first place, like $8/5$ or $5/3$.