Heron method

The Heron method , Heron's approximation method or Babylonian root extraction is a calculation method for calculating an approximation of the square root of a real number . ${\displaystyle a>0}$

Procedure

Heron method for calculating with three different starting values ${\displaystyle {\sqrt {2}}}$

The iteration equation of the Heron method can be derived from the Newton method for the zero of the quadratic function . With results from the recursion formula of Newton's method , the iteration : ${\displaystyle f(x)=x^{2}-a}$ ${\displaystyle f'(x)=2x}$ ${\displaystyle x_{n+1}=x_{n}-{\frac {f(x_{n})}{f'(x_{n})}}}$

${\displaystyle x_{n+1}=x_{n}-{\frac {x_{n}^{2}-a}{2x_{n}}}={\frac {x_{n}^{2}+a}{2x_{n}}}={\frac {1}{2}}\cdot \left(x_{n}+{\frac {a}{x_{n}}}\right)}$.

As long as it is not equal to zero, the start value of the iteration can be set arbitrarily; the iteration always converges. It should be noted that with negative starting values, the iteration converges to the negative square root. A qualified estimate for the starting value is obtained from the Taylor series expansion of the binomial series by 1, the first two terms of which provide this equation: ${\displaystyle x_{0}}$ ${\displaystyle x_{0}={\tfrac {a+1}{2}}}$

The Heron method is one of the fixed point methods . If one sets , then applies to the fixed point (which fulfills the condition ) with the (positive) solution . ${\textstyle \varphi (x)={\frac {1}{2}}\cdot \left(x+{\frac {a}{x}}\right)}$ ${\textstyle \varphi (x)=x}$ ${\textstyle x^{2}=a}$ ${\textstyle x={\sqrt {a}}}$

example

The sequence members of the Babylonian root sequence with the starting value . ${\displaystyle a_{0}={\tfrac {1}{4}}}$

In the following simple example, the square root of 9 is shown as an approximation with three calculation steps to the true value . With the start value for the iteration is calculated and inserted into the iteration rule: ${\displaystyle \textstyle {\sqrt {9}}=3}$ ${\displaystyle \textstyle a=9}$ ${\displaystyle \textstyle x_{0}={\frac {9+1}{2}}=5}$

${\displaystyle x_{1}={\frac {1}{2}}\cdot \left(5+{\frac {9}{5}}\right)={\frac {1}{2}}\cdot {\frac {34}{5}}={\frac {34}{10}}=3{,}4}$

${\displaystyle x_{2}={\frac {1}{2}}\cdot \left({\frac {34}{10}}+{\frac {9}{\frac {34}{10}}}\right)={\frac {1}{2}}\cdot \left({\frac {34}{10}}+{\frac {90}{34}}\right)={\frac {257}{85}}=3{,}0235294\dots }$

${\displaystyle x_{3}={\frac {1}{2}}\cdot \left({\frac {257}{85}}+{\frac {9}{\frac {257}{85}}}\right)={\frac {1}{2}}\cdot \left({\frac {257}{85}}+{\frac {765}{257}}\right)={\frac {65537}{21845}}=3{,}000091554\dots .}$

convergence

The procedure can be expressed as a recursively defined sequence as follows:

${\displaystyle x_{n+1}={\frac {1}{2}}\cdot \left(x_{n}+{\frac {a}{x_{n}}}\right)}$.

It is a purely positive consequence. One can now show that the -th term is for all of them . To do this, one shows the equivalent inequality : ${\displaystyle n\geq 1}$ ${\displaystyle n}$ ${\displaystyle x_{n}\geq {\sqrt {a}}}$ ${\displaystyle x_{n}^{2}-a\geq 0}$

${\displaystyle x_{n}^{2}-a={\frac {1}{4}}\cdot \left(x_{n-1}+{\frac {a}{x_{n-1}}}\right)^{2}-a={\frac {1}{4}}\cdot \left(x_{n-1}-{\frac {a}{x_{n-1}}}\right)^{2}\geq 0}$.

We also show that a monotonically decreasing sequence is: ${\displaystyle \left(x_{n}\right)}$

${\displaystyle x_{n+1}-x_{n}={\frac {1}{2}}\cdot \left(x_{n}+{\frac {a}{x_{n}}}\right)-x_{n}={\frac {a}{2x_{n}}}-{\frac {x_{n}}{2}}={\frac {a-x_{n}^{2}}{2x_{n}}}\leq 0}$.

Due to the limitedness and monotony shown, the sequence must converge, from above to the root sought:

${\displaystyle x={\frac {1}{2}}\cdot \left(x+{\frac {a}{x}}\right)\Leftrightarrow x^{2}=a\Leftrightarrow x={\sqrt {a}}}$.

Since the Heron method can be derived from Newton's approximation method and the zero to be calculated is simple, the order of convergence is 2.

The method converges very quickly if there is already a good approximation. The number of correct digits is roughly doubled with each step. However, if the first approximation is bad, it takes many steps to get a good approximation.

For example, if the root is to be calculated from an integer with 200 binary digits and you start with the first approximation, then the approximation becomes about one binary digit shorter with each step, i.e. approx. H. only after about 100 steps does the approximation have the correct length of 100 places. Then six to seven more steps ( ) are sufficient to correctly calculate all 100 places before the decimal point. ${\displaystyle a}$ ${\displaystyle x_{0}=a}$ ${\displaystyle \log _{2}(100)}$

It is therefore advisable to determine a starting value that is as precise as possible . In the example, you should first determine the bit length of and use a start value with half the length. ${\displaystyle x_{0}}$ ${\displaystyle \lfloor \log _{2}(a)\rfloor +1}$ ${\displaystyle a}$

Error estimation

The following applies to the Heron sequence : ${\displaystyle (x_{n})_{n\geq 1}}$

${\displaystyle {\frac {a}{x_{n}}}\leq {\sqrt {a}}\leq x_{n}}$ (Inclusion),

and the following estimate for the error

${\displaystyle x_{n}-{\sqrt {a}}={\frac {1}{2x_{n-1}}}\left(x_{n-1}-{\sqrt {a}}\right)^{2}\leq {\frac {1}{2{\sqrt {a}}}}\left(x_{n-1}-{\sqrt {a}}\right)^{2}}$( quadratic convergence ).

This error estimation has the disadvantage that it is not known but should be calculated. Using the above inclusion, the following workable estimate is obtained: ${\displaystyle {\sqrt {a}}}$

${\displaystyle x_{n}-{\sqrt {a}}={\frac {1}{2x_{n-1}}}\left(x_{n-1}-{\sqrt {a}}\right)^{2}\leq {\frac {1}{2x_{n-1}}}\left(x_{n-1}-{\frac {a}{x_{n}}}\right)^{2}={\frac {1}{2x_{n-1}\cdot x_{n}^{2}}}\left(x_{n-1}\cdot x_{n}-a\right)^{2}}$.

Applied to the above example one obtains:

${\displaystyle x_{3}-3={\frac {1}{2x_{2}}}\left(x_{2}-3\right)^{2}=0{,}000091554\dots \leq {\frac {1}{2x_{2}\cdot x_{3}^{2}}}\left(x_{2}\cdot x_{3}-9\right)^{2}=0{,}0000922\dots }$.

For the relative error

${\displaystyle \varepsilon _{n}={\frac {x_{n}-{\sqrt {a}}}{\sqrt {a}}}}$

the recursion applies

${\displaystyle \varepsilon _{n+1}={\frac {\varepsilon _{n}^{2}}{2(1+\varepsilon _{n})}}}$.

The consequence of is therefore independent of the starting approximation for a given relative error . ${\displaystyle \varepsilon _{n}}$ ${\displaystyle \varepsilon _{0}}$ ${\displaystyle a}$

Geometric illustration of the Heron process

The Heron method is based on the idea that a square with area has a side length of . The starting point of the procedure is any rectangle with area . Step by step, the aspect ratio of the rectangle is changed so that its shape more and more approximates that of a square, while the area remains the same. The side lengths of the rectangle are the approximate values ​​for . ${\displaystyle A}$ ${\displaystyle {\sqrt {A}}}$ ${\displaystyle A}$ ${\displaystyle {\sqrt {A}}}$

In the first step any side length is chosen for the rectangle. So that this has the desired area, the second side length is calculated using the formula ${\displaystyle x_{0}}$

${\displaystyle y_{0}={\frac {A}{x_{0}}}}$

calculated. As an example, the square root of 9 should be calculated. The value 9 is selected for one side length, so that the other side length is calculated as 1. The first rectangle therefore has the following shape.

The resemblance of this rectangle to a square is little. This is also expressed by the fact that the side lengths 1 and 9 are very poor approximations for the root of 9.

To get a better approximation of a square, the long side needs to be shortened and the short side lengthened. The new length of the long side becomes the mean

${\displaystyle x_{1}={\frac {x_{0}+y_{0}}{2}}}$

of the two previous side lengths taken. The length of the other side is calculated as above

${\displaystyle y_{1}={\frac {A}{x_{1}}}}$

In the example, the average side length is 5. The corresponding short side has a length of 1.8.

Here, too, the resemblance to a square is still small. However, the new rectangle is more compact compared to the previous one.

The process described is repeated in every further step of the Heron process. The mean value of the side lengths of a rectangle corresponds to the length of the long side of the new rectangle and the length of the short side can be calculated from this as described above. In the example, the following two rectangles are created in the next two steps.

The last rectangle is almost square. The side length of 3.024 is accordingly close to 3, the exact value of . ${\displaystyle {\sqrt {9}}}$

Implementation in software

The method is particularly suitable for implementation in software, since only basic arithmetic operations are required, so it is rarely needed today in view of the wide availability of numerical processor hardware.

If a floating point representation with a two-way exponent is also used, the approach becomes relatively simple, as an example , the root of 5 is considered and the relative error to the final value (i.e. abs ((x _i  - x) / x)) is followed:

• First an even number is split off from this two-way exponent, so that either a 0 or 1 remains as the exponent, ie the number is normalized to the interval [½, 2] . In this interval the square root function is only a slightly curved curve, so it can be treated numerically well. Example: For the  time being, only a = 1.25 with the goal x = 1.118034 will be treated. ${\displaystyle {\sqrt {5}}={\sqrt {4\cdot 1{,}25}}=2\cdot {\sqrt {1{,}25}}\approx 2\cdot 1{,}118034=2{,}236068}$
• As a starting value for the actual iteration, this curve is approximated by an even simpler one that can be calculated directly (without iteration). This initial calculation is used to determine the starting value with which the following iteration is started. This curve can be applied more or less laboriously, with the increasingly complicated approaches below, one iteration step can be saved if necessary:
  • a simple constant (for example 1), for
    example:  x ₀  = 1; rel. Error = 1.1 * 10 ⁻¹ ;
  • a straight line with a gradient of 1/2 and an additive constant of 1/2 (to simplify the following case),
    example:  x ₀ = 1/2 + 1.25 / 2 = 1.125; rel. Error = 6.2 * 10 ⁻³ ;
  • a straight line with a gradient of 1/2 and an additive, optimized constant of , example:  x ₀ = 0.929683 / 2 + 1.25 / 2 = 1.089841; rel. Error = 2.5 * 10 ⁻² ; ${\displaystyle \left(2\cdot {\sqrt[{4}]{2}}-{\sqrt {2}}\right)^{2}/2\approx 0{,}4648415}$
    
  • a straight line with an optimized slope and an additive constant (not considered here).
• A fixed number of iteration steps is carried out on the basis of the start value x ₀ determined in this way . The number required to achieve the desired accuracy can be calculated directly as the worst case within the start interval thanks to the above error estimate. For example, with a 32-bit mantissa and the middle starting point, you only need three steps. This fixed number saves considerably more time-consuming queries to achieve accuracy. Replacing the mentioned constants with the number 1.0 does not change this. Even the even more complicated approach would not save a further iteration step, at least with this accuracy. It can look different with higher accuracy requirements.
  Example with three steps according to approach 1 (constant 1, with the other approaches it converges one step faster):
  x ₁ = (x ₀ + a / x ₀ ) / 2 = (1 + 1.25 / 1) / 2 = 1.125; rel. Error = 6.2 * 10 ⁻³
  x ₂ = (x ₁ + a / x ₁ ) / 2 = (1.125 + 1.25 / 1.125) / 2 = 1.118056; rel. Error = 2.0 * 10 ⁻⁵
  x ₃ = (x ₂ + a / x ₂ ) / 2 = (1.118056 + 1.25 / 1.118056) / 2 = 1.118034; rel. Error = 0
  You can see the effect of the quadratic convergence, that the error is squared from step to step or the number of valid digits or the negative error exponent roughly doubles.
• Finally, the exponent is restored by adding half of the value split off in the first step.
  Example:  x = 2 * x ₃  = 2.236068.

Generalization of the procedure

This procedure can be generalized so that for is calculated. The larger is, the more steps are needed to accurately calculate the root. ${\displaystyle {\sqrt[{k}]{a}}}$ ${\displaystyle a>0}$ ${\displaystyle k}$

The Newton method is used to determine the positive zero of the function . The iteration rule follows from the recursion formula of Newton's method : ${\displaystyle {\sqrt[{k}]{a}}}$ ${\displaystyle f(x)=x^{k}-a}$ ${\displaystyle f'(x)=kx^{k-1}}$ ${\displaystyle x_{n+1}=x_{n}-{\frac {f(x_{n})}{f'(x_{n})}}}$

${\displaystyle x_{n+1}=x_{n}-{\frac {x_{n}^{k}-a}{kx_{n}^{k-1}}}={\frac {(k-1)x_{n}^{k}+a}{kx_{n}^{k-1}}}.}$

Here the sequence has to be started with a suitable starting value for the searched value of . ${\displaystyle x_{0}}$ ${\displaystyle {\sqrt[{k}]{a}}}$

The same convergence statements apply to positive integers as above for ${\displaystyle k}$ ${\displaystyle k=2.}$

Determination of the reciprocal

For one obtains a procedure with which (without using the division!) The reciprocal value can be calculated approximately: ${\displaystyle k=-1}$ ${\displaystyle {\sqrt[{-1}]{a}}=1/a}$

${\displaystyle x_{n+1}={\frac {(-1-1)x_{n}^{-1}+a}{(-1)x_{n}^{-1-1}}}=2x_{n}-ax_{n}^{2}=(2-ax_{n})\cdot x_{n}.}$

This procedure converges to the quadratic for all ${\displaystyle x_{0}\in \left(0,2/a\right)}$ ${\displaystyle 1/a.}$

For the first computers without built-in division, iteration made it possible to reduce this operation to multiplication and subtraction. The division of two numbers was carried out in such a way that the reciprocal of the denominator was determined and multiplied by the numerator.

example

It should be calculated approximately with the start value : ${\displaystyle 1/3}$ ${\displaystyle x_{0}={\frac {1}{2}}<{\frac {2}{3}}}$

${\displaystyle x_{1}=\left(2-3\cdot {\frac {1}{2}}\right)\cdot {\frac {1}{2}}={\frac {1}{4}}=0{,}25,}$

${\displaystyle x_{2}=\left(2-3\cdot {\frac {1}{4}}\right)\cdot {\frac {1}{4}}={\frac {5}{16}}=0{,}3125,}$

${\displaystyle x_{3}=\left(2-3\cdot {\frac {5}{16}}\right)\cdot {\frac {5}{16}}={\frac {85}{256}}=0{,}33203125.}$

For the start value one gets ${\displaystyle x_{0}={\frac {2}{3}}}$

${\displaystyle x_{1}=\left(2-3\cdot {\frac {2}{3}}\right)\cdot {\frac {2}{3}}=0,}$

${\displaystyle x_{2}=\left(2-3\cdot 0\right)\cdot 0=0,}$

thus no convergence to the sought value of ${\displaystyle {\frac {1}{3}}.}$

Historical

The method was already known in Mesopotamia at the time of Hammurapi I (approx. 1750 BC), a king of Babylon . It was described by Heron of Alexandria in the first book of his work Metrica around 100 AD .

literature

• Bernd Ziegenbalg: Algorithms: From Hammurapi to Gödel . Harri Deutsch Verlag 2007, ISBN 978-3-8171-1814-4 , pp. 54-59.
• David Fowler, Eleanor Robson : Square Root Approximations in Old Babylonian Mathematics (PDF; 215 kB). Historica Mathematica 25, 1998, pp. 366-378.
• Hans R. Schwarz, Norbert Köckler: Numerical Mathematics. 6th edition. Teubner, Stuttgart 2006, ISBN 3-519-42960-8 , pp. 189-192.

Web links

• The Heron method for calculating the roots on arndt-bruenner.de (explanation and example calculator)

Remarks

1. ↑ Start value: If the output value is already available as a binary number (in the place value system), you can simply count where its highest value '1' is; Start value is then . If the output value is in (binary) exponential notation, the exponent can simply be halved as the start value (shift by 1 bit to the right). See also section Implementation in software ${\displaystyle i}$ ${\displaystyle 1\cdot 2^{i/2}}$

Individual evidence

1. ↑ Suitable transformations: zeros and fixed points. In: Montanuniversität Leoben . February 23, 2005, accessed August 27, 2019 . 
2. ↑ Bernd Ziegenbalg: Algorithms: From Hammurapi to Gödel . Harri Deutsch Verlag 2007, ISBN 978-3-8171-1814-4 , p. 54 ( excerpt (Google) )
3. ↑ John J. O'Connor, Edmund F. Robertson :  Heron method. In: MacTutor History of Mathematics archive .