Methods of computing square roots

This article presents and explains several methods which can be used to calculate square roots.

Exponential identity

Pocket calculators typically implement good routines to compute the exponential function and the natural logarithm, and then compute the square root of r using the identity

${\sqrt {r}}=e^{{\frac {1}{2}}\ln r}$.

The same identity is used when computing square roots with logarithm tables or slide rules.

Rough estimation

Many of the methods for calculating square roots require an initial seed value. If the initial value is too far from the actual square root, the calculation will be slowed down. It is therefore useful to have a rough estimate, which may be very inaccurate but easy to calculate. One way of obtaining such an estimate for ${\sqrt {r}}$ is calculating 3^D, where D is the number of digits in r (if r < 1, D will be negative).

A better method of rough estimation is this:

If D is odd, D = 2n + 1, then use ${\sqrt {r}}\approx 2\cdot 10^{n}$

If D is even, D = 2n + 2, then use ${\sqrt {r}}\approx 6\cdot 10^{n}$

When working in the binary numeral system (as computers do internally), an alternative method is to use $2^{\left\lfloor D/2\right\rfloor }$ (here D is the number of binary digits).

Babylonian method

A commonly used algorithm for approximating ${\sqrt {r}}$ is known as the "Babylonian method" and can be derived from (but predates) Newton's method. This is a quadratically convergent algorithm, which means that the number of correct digits of $r$ roughly doubles with each step. It proceeds as follows:

1. Start with an arbitrary positive start value $x$ (the closer to the root, the better).
2. Replace $x$ by the average of $x$ and r/x.
3. Repeat steps 2 and 3, until the desired accuracy is achieved.

It can also be represented as:

$x_{0}\approx {\sqrt {r}}$

$x_{n+1}={\frac {1}{2}}\left(x_{n}+{\frac {r}{x_{n}}}\right)$

$\lim _{n\to \infty }x_{n}={\sqrt {r}}$.

This algorithm works equally well in the p-adic numbers, but cannot be used to identify real square roots with p-adic square roots; it is easy, for example, to construct a sequence of rational numbers by this method which converges to $+3$ in the reals, but to $-3$ in the 2-adics.

Example

We'll calculate ${\sqrt {r}}$, where r = 125348, to 6 significant figures. Here D, the number of digits in r, is 6. Then:

$x_{0}=3^{D}=3^{6}=729.000\,\!$

$x_{1}={\frac {1}{2}}\left(x_{0}+{\frac {r}{x_{0}}}\right)={\frac {1}{2}}\left(729.000+{\frac {125348}{729.000}}\right)=450.472$

$x_{2}={\frac {1}{2}}\left(x_{1}+{\frac {r}{x_{1}}}\right)={\frac {1}{2}}\left(450.472+{\frac {125348}{450.472}}\right)=364.365$

$x_{3}={\frac {1}{2}}\left(x_{2}+{\frac {r}{x_{2}}}\right)={\frac {1}{2}}\left(364.365+{\frac {125348}{364.365}}\right)=354.191$

$x_{4}={\frac {1}{2}}\left(x_{3}+{\frac {r}{x_{3}}}\right)={\frac {1}{2}}\left(354.191+{\frac {125348}{354.191}}\right)=354.045$

$x_{5}={\frac {1}{2}}\left(x_{4}+{\frac {r}{x_{4}}}\right)={\frac {1}{2}}\left(354.045+{\frac {125348}{354.045}}\right)=354.045$

Therefore, ${\sqrt {125348}}\approx 354.045$.

Bakhshali approximation

This is a method for finding an approximation to a square root which was described in an ancient manuscript known as the Bakhshali Manuscript. It is equivalent to the Babylonian method applied twice. The original presentation goes as follows: To calculate ${\sqrt {r}}$, let N² be the nearest perfect square to r. Then, calculate:

$d=r-N^{2}\,\!$

$P={\frac {d}{2N}}$

$A=N+P\,\!$

${\sqrt {r}}\approx A-{\frac {P^{2}}{2A}}$

This can be also written as:

${\sqrt {r}}\approx {\frac {N^{4}+6N^{2}r+r^{2}}{4N^{3}+4Nr}}$

${\sqrt {N^{2}+d}}\approx N+{\frac {d}{2N}}-{\frac {d^{2}}{8N^{3}+4Nd}}$

Example

We'll find ${\sqrt {9.2345}}$.

$N=3\,\!$

$d=9.2345-3^{2}=0.2345\,\!$

$P={\frac {0.2345}{2\cdot 3}}=0.0391$

$A=3+0.0391=3.0391\,\!$

${\sqrt {9.2345}}\approx 3.0391-{\frac {0.0391^{2}}{2\cdot 3.0391}}=3.0388$

Digit by digit calculation

This is a method to find each digit of the square root in a sequence. It is much slower than the Babylonian method, but it has several advantages:

• It can be easier for manual calculations.
• Every digit found is known to be correct.
• If the square root has an expansion which terminates, the algorithm will terminate after the last digit is found. It can thus be used to check whether a given integer is a square number.

Napier's bones include an aid for the execution of this algorithm. The Shifting nth-root algorithm is a generalization of this method.

The algorithm works for any base, and naturally, the way it proceeds depends on the base chosen. We will describe the method for the decimal system. It goes as follows:

Write the the decimal expansion of the number. The numbers are laid out in a fashion similar to the long division algorithm, and the final result will appear above the original number. Divide the digits of the number into pairs, starting from the decimal point. Each pair of digits corresponds to a single digit of the result.

For each iteration:

1. Bring down the most significant (leftmost) pair of digits not yet used (if all the digits of the number have been used, use the digits 00) and append them to the remainder of the previous step (for the first step, the remainder is 0). This will be the "current value".
2. Let s denote the part of the result found so far, ignoring any decimal point (for the first step, s = 0), determine the greatest digit x such that $y=(20s+x)x$ does not exceed the current value (20s + x is simply twice s, with the digit x appended to the right). Place the digit $x$ as the next digit of the result.
3. Subtract $y$ from the current value to form a new remainder.
4. If the remainder is zero and there are no more digits to bring down the algorithm has terminated. Otherwise continue with step 1.

Examples

Find the the square root of 152.2756.

          1  2. 3  4 
     \/  01 52.27 56                            

x        01                   1*1 <= 1 < 2*2                 x = 1
         01                     y = x*x = 1*1 = 1
2x       00 52                22*2 <= 52 < 23*3              x = 2
         00 44                  y = 2x*x = 22*2 = 44                      
24x         08 27             243*3 <= 827 < 244*4           x = 3       
            07 29               y = 24x*x = 243*3 = 729
246x           98 56          2464*4 <= 9856 < 2465*5        x = 4       
               98 56            y = 246x*x = 2464*4 = 9856
               00 00          Algorithm terminates: Answer is 12.34

Find the square root of 2.

          1. 4  1  4  2
     \/  02.00 00 00 00

x        02                   1*1 <= 2 < 2*2                 x = 1
         01                     y = x*x = 1*1 = 1
2x       01 00                24*4 <= 100 < 25*5             x = 4
         00 96                  y = 2x*x = 24*4 = 96                      
28x         04 00             281*1 <= 400 < 282*2           x = 1       
            02 81               y = 28x*x = 281*1 = 281
282x        01 19 00          2824*4 <= 11900 < 2825*5       x = 4       
            01 12 96            y = 282x*x = 2824*4 = 11296
2828x          06 04 00       28282*2 <= 60400 < 28283*3     x = 2
                              Desired accuracy achieved: Answer is 1.4142

Iterative methods for inverse square roots

The following are iterative methods for finding the inverse square root of r, ${\frac {1}{\sqrt {r}}}$. Once it has been found, we can find ${\sqrt {r}}$ by simple multiplication: ${\sqrt {r}}={\frac {1}{\sqrt {r}}}\cdot r$. The iterations involve only multiplication, and not division; They are therefore faster than the Babylonian method. However, they are not stable; If the initial value is not close to the root, the iterations will diverge away from it rather than converge to it.

• This is found by applying Newton's method to the equation ${\frac {1}{x^{2}}}-r=0$:

$x_{n+1}={\frac {1}{2}}x_{n}(3-rx_{n}^{2})$

• This requires less iterations to converge, but involves more operations per iteration:

$y_{n}=rx_{n}^{2}\,\!$

$x_{n+1}={\frac {1}{8}}x_{n}(15-y_{n}(10-3y_{n}))$

Taylor series

If N is an approximation to ${\sqrt {r}}$, a better approximation can be found by using the Taylor series of the square root function:

${\sqrt {N^{2}+d}}=\sum _{n=0}^{\infty }{\frac {(-1)^{n}(2n)!d^{n}}{(1-2n)n!^{2}4^{n}N^{2n-1}}}=N+{\frac {d}{2N^{2}}}-{\frac {d^{2}}{8N^{3}}}+{\frac {d^{3}}{16N^{5}}}-{\frac {5d^{3}}{128N^{7}}}+\dots$

As an iterative method, the order of convergence is equal to the number of terms used. With 2 terms, it is identical to the Babylonian method; With 3 terms, each iteration takes almost as many operations as the Bakhshali approximation, but converges more slowly. Therefore, this is not a particularly efficient way of calculation.

Continued fraction expansion

Quadratic irrationals (numbers of the form ${\frac {a+{\sqrt {b}}}{c}}$, where a, b and c are integers), and in particular, square roots of integers, have periodic continued fractions. Sometimes we may be interested not in finding the numerical value of a square root, but rather in its continued fraction expansion. The following iterative algorithm can be used for this purpose (r is any natural number which is not a perfect square):

$m_{0}=0\,\!$

$d_{0}=1\,\!$

$a_{0}=\left\lfloor {\sqrt {r}}\right\rfloor \,\!$

$m_{n+1}=d_{n}a_{n}-m_{n}\,\!$

$d_{n+1}={\frac {r-m_{n+1}^{2}}{d_{n+1}}}\,\!$

$a_{n+1}=\left\lfloor {\frac {a_{0}+m_{n+1}}{d_{n+1}}}\right\rfloor \,\!$

The algorithm terminates when the triplet (m_n, d_n, a_n) is identical to one encountered before; The expansion will repeat from then on. The sequence [a₀; a₁, a₂, a₃, …] is the continued fraction expansion:

${\sqrt {r}}=a_{0}+{\cfrac {1}{a_{1}+{\cfrac {1}{a_{2}+{\cfrac {1}{a_{3}+\,\cdots }}}}}}$

Pell's equation

Pell's equation and its variants yield a method for efficiently finding continued fraction convergents of square roots of integers. However, it can be complicated to execute, and usually not every convergent is generated. The ideas behind the method are as follows:

• If (p, q) is a solution (where p and q are integers) to the equation $p^{2}=rq^{2}\pm 1$, then ${\frac {p}{q}}$ is a continued fraction convergent of ${\sqrt {r}}$, and as such, is an excellent rational approximation to it.
• If (p_a, q_a) and (p_b, q_b) are solutions, then so is:

$p=p_{a}p_{b}+rq_{a}q_{b}\,\!$

$q=p_{a}q_{b}+p_{b}q_{a}\,\!$

• More generally, if (p₁, q₁) is a solution, then it is possible to generate a sequence of solutions (p_n, q_n) satisfying:

$p_{m+n}=p_{m}p_{n}+rq_{m}q_{n}\,\!$

$q_{m+n}=p_{m}q_{n}+p_{n}q_{m}\,\!$

The method is as follows:

• Find natural numbers (0 excluded) p₁ and q₁ such that $p_{1}^{2}=rq_{1}^{2}\pm 1$. This is the hard part; It can be done either by guessing, or by using fairly sophisticated techniques.

• To generate a long list of convergents, iterate:

$p_{n+1}=p_{1}p_{n}+rq_{1}q_{n}\,\!$

$q_{n+1}=p_{1}q_{n}+p_{n}q_{1}\,\!$

• To find the larger convergents quickly, iterate:

$p_{2n}=p_{n}^{2}+rq_{n}^{2}\,\!$

$q_{2n}=2p_{n}q_{n}\,\!$

• In either case, ${\frac {p_{n}}{q_{n}}}$ is a rational approximation satisfying:

$\left|{\frac {p_{n}}{q_{n}}}-{\sqrt {r}}\right|<{\frac {1}{q_{n}^{2}{\sqrt {r}}}}$

Approximations that depend on IEEE representation

On computers, a very rapid Newton's method based approximation to the square root can be obtained for floating point numbers when computers use an IEEE (or sufficiently similar) representation.

float fastsqrt(float val) {
    int tmp = *(int *)&val;
    tmp -= 1<<23; /* Remove IEEE bias from exponent (-2^23) */
    /* tmp is now an appoximation to logbase2(val) */
    tmp = tmp >> 1; /* divide by 2 */
    tmp += 1<<23; /* restore the IEEE bias from the exponent (+2^23) */
    return *(float *)&tmp;
}

In the above, the operations to add and remove the IEEE bias can be combined into a single operation. An additional adjustment can be made in the same operation to reduce the maximum relative error. So, the three operations, not including the cast, can be rewritten as:

   tmp = (1<<23) + (tmp >> 1) - 1<<22 + m;

Where m is some magic number that corresponds to a calculation involving an initial guess used in the Newton's approximation.

A variant of the above routine can be used to compute the inverse square root $x^{-{1 \over 2}}$ instead, which has been attributed to the programmers of the game Doom. As an approximation, it produced a relative error of less than 4%, and in computer graphics it is a very efficient way to normalize a vector.

See also

• Alpha max plus beta min algorithm
• Integer square root
• Shifting nth-root algorithm
• N-th root algorithm

External links

• Weisstein, Eric W. "Square root algorithms". MathWorld.