Methods of computing square roots

There are many methods to compute square roots. Below is a list and explanation of methods that have been discovered throughout the history of human mathematical study.

Exponential identity

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

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

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

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. 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.

This is a quadratically convergent algorithm, which means that the number of correct digits of $r$ roughly doubles with each step. It can also be represented as:

$x_{n+1}=0.5(x_{n}+{\frac {r}{x_{n}}})$

$\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.

Rough estimation

Many of the methods for calculating square roots , such as the Babylonian method, 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).

Bakhshali approximation

The Bakhshali square root approximation is a mathematical method for finding an approximation to a square root which was described in an ancient manuscript, written in ancient India sometime between the 200 BC and the 400 CE, known as the Bakhshali Manuscript because it was discovered in 1881 near the village of Bakhshali (or Bakhshalai) in the Yusufzai subdivision of the Peshawar district (now part of Pakistan). The manuscript was made out of the leaves of birch bark tree and was written in the Sarada script.

The Bakhshali square root approximation is just the simple approximation applied twice. This gives the same result as using the simple approximation as an initial value to Newton's method, then doing just one step of the Newton's method. It is just a little slower to compute than Newton's.

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

Let $A=N+P$

${\sqrt {N^{2}+d}}\approx A-{\frac {P^{2}}{2A}}$

When expanded, the equation becomes

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

Example

Find ${\sqrt {9.2345}}$

Using $N=3$ and $d=9.2345-3^{2}=0.2345$

${\sqrt {3^{2}+d}}\approx 3+{\frac {d}{6}}-{\frac {d^{2}}{216+12\,d}}$

${\sqrt {3^{2}+0.2345}}\approx 3+{\frac {0.2345}{6}}-{\frac {0.2345^{2}}{216+12\cdot 0.2345}}$

${\sqrt {3^{2}+0.2345}}\approx 3+0.039083-{\frac {0.055}{216+2.814}}$

${\sqrt {3^{2}+0.2345}}\approx 3+0.039083-0.000251$

${\sqrt {3^{2}+0.2345}}\approx 3.038832$

We can get to the same result faster if we do not use the expansion:

$P=0.2345/(2*3)=0.039083$

$A=3+P=3.039083$

$result=3-0.039083^{2}/(2*A)=3.038832$

And we can save one multiplication if instead of using this last line we compute the result directly by Newton's formula, like this:

$result=(A+9.2345/A)/2=3.038832$

this always gives the exact same result.

An exact long-division-like algorithm

This method, while much slower than the Babylonian method, has the advantage that it is exact: if the given number has a square root whose decimal representation terminates, then the algorithm terminates and produces the correct square root after finitely many steps. It can thus be used to check whether a given integer is a square number.

Write the number in decimal and divide it into pairs of digits starting from the decimal point. The numbers are laid out similar to the long division algorithm and the final square root will appear above the original number.

For each iteration:

1. Bring down the most significant pair of digits not yet used and append them to any remainder. This is the current value referred to in steps 2 and 3.
2. If $s$ denotes the part of the result found so far, determine the greatest digit $x$ that does not make $y=x(20s+x)$ exceed the current value. Place the new digit $x$ on the quotient line.
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.

Example: What is the square root of 152.2756?

           1  2. 3  4 
       |  01 52.27 56                              1    The digit 1 is a solution digit
x         01                   1(20*0+1)=1        +1 
          00 52                                    22   The digit 2 is a solution digit
2x        00 44                2(20*1+2)=44       + 2 
             08 27                                 243  The digit 3 is a solution digit
24x          07 29             3(20*12+3)=729     +  3 
                98 56                              2464 The digit 4 is a solution digit
246x            98 56          4(20*123+4)=9856       4
                00 00          Algorithm terminates: answer is 12.34

Although demonstrated here for base 10 numbers, the procedure works for any base, including base 2. In the description above, 20 means double the number base used, in the case of binary this would really be 100. The algorithm is in fact much easier to perform in base 2, as in every step only the two digits 0 and 1 have to be tested. Napier's bones used this algorithm. See also Shifting nth-root algorithm.

More accurate rough estimation

Similar to the rough estimation procedure above. This version is slightly more accurate because it also uses the additional information provided by the first digit of the number r. The steps are almost the same as the one for previous rough estimate. This estimate can then be improved by any iterative method.

Provided that r ≫ 1

Steps

1. Take the integer part of the number r. Z = int(r)
2. Count the number of digits in Z. Let D be the number of digits.
3. Calculate the value of 3.16^D. Note: 3.16 ≈ √10
4. The estimate is the value obtained in step 3 multiplied by an adjustment factor based on the first (left most) digit of the number n:
   E = ADJ * 3.16^D

Table of ADJ

First Digit	1	2	3	4	5	6	7	8	9
ADJ	0.382	0.497	0.590	0.670	0.741	0.806	0.866	0.922	0.974

where the derivation of the adjustment is the average of the square root of first digit and the square root of the first digit incremented by one, all divided by √10.

The accuracy of this method is very dependant on the number's first digit. Due to logrithmic spacings, errors for numbers starting with a "1" can range up to 18 percent, while errors for numbers starting with a "9" are less than 3 percent.

Example

r	Z	D	3.16 ^ D	ADJ	Estimate (E)	Actual Value of √r
723.47	723	3	31.55	0.866	27.32	26.90
5396.37	5396	4	99.71	0.741	73.89	73.46
24956.41	24956	5	315.09	0.497	156.60	157.98
789345.464	789345	6	995.7	0.866	862.28	888.45
198764.389	198764	6	995.7	0.382	380.36	445.83

Inverse square root approximation

There is an algorithm for finding the inverse square root $\left({\frac {1}{\sqrt {r}}}\right)$ very quickly.

Because ${\sqrt {r}}=r\times \left({\frac {1}{\sqrt {r}}}\right)$, if we know the value of the inverse square root, we can easily get the value for the square root.

Steps for inverse square root approximation

Step 1 Start the initial value $X_{0}$ an approximate value of the inverse square root which is accurate to 2 significant digits.

$X_{0}=\left({\frac {1}{\sqrt {r}}}\right)_{\mbox{approx}}$

Step 2 Calculate the next iteration with the following algorithm.

$X_{n+1}=X_{n}\,\left({\frac {15}{8}}-Y_{n}\,\left({\frac {5}{4}}-{\frac {3}{8}}\,Y_{n}\right)\right)$ where $Y_{n}=r\,X_{n}^{2}$

Step 3 Finally obtain the value of ${\sqrt {r}}$ with

${\sqrt {r}}=r\,X_{infinity}$

Because $X_{infinity}$ would converge to the value of the inverse square root of r.

Example

Find the ${\sqrt {7}}$

Step 1 Find initial value $X_{0}$ using Bakhshali square root approximation.

${\sqrt {7}}={\sqrt {2^{2}+3}}={\sqrt {N^{2}+d}}$

$P\,=\,{\frac {d}{2N}}\,=\,{\frac {3}{4}}$

$A=N+P\,=\,{\frac {11}{4}}$

${\sqrt {N^{2}+d}}\approx A-{\frac {P^{2}}{2A}}\approx {\frac {12500}{4721}}\approx 2.64773$

$X_{0}={\frac {1}{2.64773}}=0.37768$

Step 2 Calculate the next iteration with the following algorithm.

$X_{n+1}=X_{n}\,\left({\frac {15}{8}}-Y_{n}\,\left({\frac {5}{4}}-{\frac {3}{8}}\,Y_{n}\right)\right)$ where $Y_{n}=r\,X_{n}^{2}$

$Y_{0}=7*X_{0}^{2}=7*0.37768^{2}=0.9984952768$

$X_{1}=0.377964472606588$

$Y_{1}=7*X_{1}^{2}=7*0.377964472606588^{2}=0.99999999786943356769463711$

$X_{2}=0.3779644730092272272145165350918641$

$Y_{2}=7*X_{2}^{2}=7*0.3779644730092272272145165350918641^{2}=0.99999999786943356769463711$

$X_{3}=0.37796447300922722721451653623418$

Step 3 Finally obtain the value of ${\sqrt {r}}$

${\sqrt {7}}\approx 7*X_{3}\approx 2.64575131106459059050161575363926$

Properties

The greatest weakness of this algorithm is that it is not stable. The initial value $X_{0}$ must be near the actual value for the algorithm to converge correctly. If the initial value is not close enough, the algorithm will diverge away from the actual value.

However when the algorithm converges, it converges quickly. 20 digits of accuracy can be obtained in 6 iterations.

Square roots using Newton iteration

Basic Newton iteration finds a single root of a function $f(x)$ given a sufficiently precise approximation to the root. The nature of which root will be given based on an approximation is dependent on the Newton fractal. The basic iteration is given by:

$x_{n+1}=x_{n}-{f(x_{n}) \over f^{\prime }(x_{n})}$.

There are two widely used functions $f(x)$ and $g(x)$ used to find the square root of a number denoted by "z".

(See also: Newton's method, Integer square root.)

First method

The first method finds the square root of "z"

$f(x)=x^{2}-z\,$

Note that both ${\sqrt {z}}$ and $-{\sqrt {z}}$ are roots of the function $f(x)$. ie $f({\sqrt {z}})=0$.

The first derivative of $f(x)$ is $f^{\prime }(x)=2x$

Thus iteration for $x_{n+1}$ is derived where:

$x_{0}=1\,\,\,$        and

$x_{n+1}\,\!$	$=x_{n}-{f(x_{n}) \over f^{\prime }(x_{n})}$
	$=x_{n}-{(x_{n}^{2}-z) \over 2x_{n}}$
	$=x_{n}-{\frac {x_{n}}{2}}+{\frac {z}{2\,\,\,x_{n}}}$
	$={\frac {x_{n}}{2}}+{\frac {z}{2\,\,\,x_{n}}}$.

So this simplifies to the Babylonian method described above.

Second method

The second method finds the reciprocal of the square root of "z".

$g(x)={\frac {1}{x^{2}}}-z$.

The two roots to $g(x)$ are ${\frac {1}{\sqrt {z}}}$ and ${\frac {-1}{\sqrt {z}}}$.

The derivative of $g(x)$ is $g^{\prime }(x)=-2x^{-3}$.

Thus iteration for $x_{n+1}$ is derived where:

$x_{0}=0.5$        and

$x_{n+1}\,\!$	$=x_{n}-{g(x_{n}) \over g^{\prime }(x_{n})}$
	$=x_{n}-{x_{n}^{-2}-z \over -2x_{n}^{-3}}$
	$=x_{n}-(-1/2){x_{n}^{3}}(x_{n}^{-2}-z)$
	$=x_{n}+(1/2)(x_{n}-zx_{n}^{3})$
	$={\frac {3x_{n}-zx_{n}^{3}}{2}}$
	$=1.5\,x_{n}-0.5\,zx_{n}^{3}=0.5\,x_{n}\,\,(3-zx_{n}^{2})$.

The square root of "z" can be found by multiplying the number just computed by "z", so this method does not require any division at all.

Example

Find the ${\sqrt {7}}$ using both methods.

$z=7\,$ Because we are looking for the square root of 7

$f(x)\,$		$g(x)\,$
$x_{0}\,$	$1\,$	$x_{0}\,$	$0.5\,$
$x_{1}\,$	${\frac {1}{2}}+{\frac {7}{2\times 1}}=4$	$x_{1}\,$	$0.5\times 0.5\,\,(3-7(0.5)^{2})=0.3125$	$x_{1}z=2.188\,$
$x_{2}\,$	${\frac {4}{2}}+{\frac {7}{2\times 4}}=2.875$	$x_{2}\,$	$0.5\times 0.3125\,\,(3-7(0.3125)^{2})=0.3619$	$x_{2}z=2.534\,$
$x_{3}\,$	${\frac {2.875}{2}}+{\frac {7}{2\times 2.875}}=2.654$	$x_{3}\,$	$0.5\times 0.3619\,\,(3-7(0.3619)^{2})=0.377$	$x_{3}z=2.639\,$
$x_{4}\,$	${\frac {2.654}{2}}+{\frac {7}{2\times 2.654}}=2.646$	$x_{4}\,$	$0.5\times 0.377\,\,(3-7(0.377)^{2})=0.378$	$x_{4}z=2.646\,$
${\sqrt {7}}\approx 2.646$		${\sqrt {7}}\approx 2.646$

Comparison

The iteration for $f(x)$ involves a division which is more time consuming than a multiplication (usually by a factor of 3 or more) in computer integer arithmetic. The iteration for $g(x)$ involves no division and is thus recommended for large integers z. However, the $f(x)$ iteration always converges, whereas the convergence of the $g(x)$ iteration is sensitive to the initial seed.

Pell's equation

Pell's equation yields a method for finding rational approximations of square roots of integers.

Finding square roots using mental arithmetic

Based on Pell's equation there is a method to calculate square roots of integers simply by subtracting odd numbers, adapted therefore for mental arithmetic.

Example: For ${\sqrt {27}}$ we start with the following sequence:

27 - 1 = 26

26 - 3 = 23

23 - 5 = 18

18 - 7 = 11

11 - 9 = 2

Five steps have been taken and thus the integer part of the square root of 27 is 5. We do not proceed any further because the sixth subtraction would give a negative answer.

Now take the rest (2) and multiply it by 100 to get the starting number for the next step. Take the answer we have already got (5) and multiply it by 20, then add 1 to get the first odd number we subtract by. This gives us 2 × 100 = 200 as the starting number and 5 × 20 + 1 = 101 as the first odd number. Subtract the successive odd numbers, 103, 105, etc. until we get to a step where the next subtraction would result in a negative number. So,

200 - 101 = 99

and 99 is less than 103 so this is the place to stop and since we only took one step we get that the next digit is 1.

For the next step the starting number will be 99 × 100 = 9900 and the answer we have already got is 51 so the first odd number will be 51 × 20 + 1 = 1021

9900 - 1021 = 8879

8879 - 1023 = 7856

7856 - 1025 = 6831

6831 - 1027 = 5804

5804 - 1029 = 4775

4775 - 1031 = 3744

3744 - 1033 = 2711

2711 - 1035 = 1676

1676 - 1037 = 639

We took nine steps so the next digit is 9.

Continuing with this procedure, the next starting number will be 639 × 100 = 63900 and the first odd number will be 519 × 20 + 1 = 10381

63900 - 10381 = 53519
53519 - 10383 = 43136
43136 - 10385 = 32751
32751 - 10387 = 22364
22364 - 10389 = 11975
11975 - 10391 = 1584

We took six steps so the next digit is 6.

The result gives us 5.196 as an approximation of the square root of 27.

Approximating with a random method for arbitrarily large integers

The idea of this method for finding the integer part of the square root of $n$ is to choose a random integer in the range $[1,n]$, and square it. If it is greater than $n$ and smaller than the right bound of the range, it becomes the new right bound. If the number squared is less than n and greater than the left bound of the range, it becomes the new left bound. In this way a computer can find a good approximation to the square root for any size number very quickly. Here is Maple code that does the job:

 intsqrt := proc(z) local A,B,a,b,c,e,r;
  A := 1;
  B:= z;
  e := z-1;
  while (e>1) do
   r := rand(A+1..B-1):
   a := r(e);
   if (a > A and a^2 < z) then
    A := a;
   end;
   if (a < B and a^2 > z) then
    B := a;
   end;
   e := B-A;
  end;
 printf("Integer part of sqrt is %d.\n", A);
 end proc:

Finding square roots using basic arithmetic

If you have access to a calculator capable of addition, subtraction, multiplication and division then the following method can be used to find the square root of a number. Note that this algorithm does not always terminate when attempting to determine the next digit. For example: 2, 3 and 7 do not terminate on the first digit.

Find ${\sqrt {r}}$

$A_{base}=r\qquad \qquad \qquad \qquad at\quad the\quad very\quad start$

$B_{base}=SOLN\times 20+1\qquad where\quad SOLN=0\quad at\quad the\quad very\quad start$

The starting values are

${\begin{matrix}D_{0}&=&9\\A_{0}&=&A_{base}+D_{0}\,(D_{0}-1)\\B_{0}&=&B_{base}+2\,(D_{0}-1)\end{matrix}}$

Now find new values for each variable

$A_{n}=A_{base}+D_{n}\,(D_{n}-1)$

$B_{n}=B_{base}+2\,(D_{n}-1)$

$D_{n+1}=floor(A_{n}\div B_{n})$

We found a digit when the value of    “$D_{n}$”    is equal to the value of    “$floor(A_{n}\div B_{n})$”

After a new digit has been found.

New $A_{base}=100\,(A_{n}-B_{n}\times D_{n})$

New $B_{base}=SOLN\times 20+1$

The step are shown in the table below

Find ${\sqrt {27}}$
$A_{base}=Z=27$ and $B_{base}=0\times 20+1=1$
$n$	$D_{n}$	$A_{n}$	$B_{n}$	$floor(A_{n}\div B_{n})$	Remarks	SOLN
n=0	$9$	$27+9\times 8=99$	$1+2\times 8=17$	$5$
n=1	$5$	$27+5\times 4=47$	$1+2\times 4=9$	$5$	STOP	$5$
$A_{base}=100\,(47-9\times 5)=200$ and $B_{base}=5\times 20+1=101$
n=0	$9$	$200+9\times 8=272$	$101+2\times 8=117$	$2$
n=1	$2$	$200+2\times 1=202$	$101+2\times 1=103$	$1$
n=2	$1$	$200+1\times 0=200$	$101+2\times 0=101$	$1$	STOP	$51$
$A_{base}=100\,(200-101\times 1)=9900$ and $B_{base}=51\times 20+1=1021$
n=0	$9$	$9900+9\times 8=9972$	$1021+2\times 8=1037$	$9$	STOP	$519$
$A_{base}=100\,(9972-1037\times 9)=63900$ and $B_{base}=519\times 20+1=10381$
n=0	$9$	$63900+9\times 8=63972$	$10381+2\times 8=10397$	$6$
n=1	$6$	$63900+6\times 5=63930$	$10381+2\times 5=10391$	$6$	STOP	$5196$
$A_{base}=100\,(63930-10391\times 6)=158400$ and $B_{base}=5196\times 20+1=103921$
n=0	$9$	$158400+9\times 8=158472$	$103921+2\times 8=103937$	$1$
n=1	$1$	$158400+1\times 0=158400$	$103921+2\times 0=103921$	$1$	STOP	$51961$
${\sqrt {27}}\approx 5.1961$

Continued fraction methods

Quadratic irrationals, that is numbers involving square roots in the form (a + √b)/c, have periodic continued fractions. This makes them easy to calculate recursively given the period. For example, to calculate √2, we make use of the fact that √2 − 1 = [0; 2, 2, 2, 2, 2, ...], and use the recurrence relation

a_{n + 1} = 1/(2 + a_n) with a₀ = 0

to obtain √2 − 1 to some specific precision specified through n levels of recurrence, and add 1 to the result to obtain √2.

Steps for finding continued fractions

Find ${\sqrt {r}}$ using continued fractions

${\sqrt {r}}=N_{0}+{\frac {1}{N_{1}+{\frac {1}{N_{2}+{\frac {1}{N_{3}+\,\cdots }}}}}}$

Let ${\sqrt {r}}={\sqrt {N_{0}^{2}+d}}=N_{0}+{\frac {1}{X_{0}}}\qquad where\quad X_{0}>1\qquad {\mbox{ and }}\quad N_{0}\in \mathbb {Z^{+}}$

thus

${\sqrt {r}}=N_{0}+{\frac {1}{X_{0}}}$

$X_{0}=N_{1}+{\frac {1}{X_{1}}}$

$X_{1}=N_{2}+{\frac {1}{X_{2}}}$

$X_{2}=N_{3}+{\frac {1}{X_{3}}}$

Step 0. We shall assume that r is not a perfect square. In other words:

${\sqrt {r}}={\sqrt {N_{0}^{2}+d}}\qquad {\mbox{ and }}\quad N_{0}\in \mathbb {Z^{+}} \quad {\mbox{ and }}0<d<1$

Step 1. Find $N_{0}$ using some other method. The best method is $N_{0}=intsqrt(r)$ using some other algorithm to determine the integer square root.

$N_{0}=intsqrt(r)$

Step 2. Find the highest lower bound (L) and lowest upper bound (U) for ${\sqrt {r}}$ where both (L) and (U) are integers.

$L\qquad <\qquad {\sqrt {r}}\qquad <\qquad U$

Hence

$L\qquad =\qquad N_{0}$

$U\qquad =\qquad N_{0}+1$

Step 3. Write $X_{0}$ in terms of ${\sqrt {r}}$.

$X_{0}={\frac {1}{{\sqrt {r}}-N_{0}}}\qquad from\qquad {\sqrt {r}}=N_{0}+{\frac {1}{X_{0}}}$

$X_{0}={\frac {1}{{\sqrt {r}}-N_{0}}}\times {\frac {{\sqrt {r}}+N_{0}}{{\sqrt {r}}+N_{0}}}={\frac {{\sqrt {r}}+N_{0}}{r-N_{0}^{2}}}={\frac {{\sqrt {r}}+N_{0}}{d}}={\frac {{\sqrt {r}}+M_{0}}{D_{0}}}$

${\mbox{ where }}\quad M_{0}=N_{0}$

and

${\mbox{ where }}\quad D_{0}=d$

$X_{0}={\frac {{\sqrt {r}}+M_{0}}{D_{0}}}$

$\downarrow$

${\mbox{ substitute }}\quad {\sqrt {r}}\quad {\mbox{ with }}\quad L\quad {\mbox{ or }}\quad U$

$\downarrow$

$X_{0LowerBound}={\frac {{\sqrt {r}}+M_{0}}{D_{0}}}={\frac {L+M_{0}}{D_{0}}}={\frac {N_{0}+M_{0}}{D_{0}}}\qquad {\mbox{ Calc the numeric value }}$

and

$X_{0UpperBound}={\frac {{\sqrt {r}}+M_{0}}{D_{0}}}={\frac {U+M_{0}}{D_{0}}}={\frac {(N_{0}+1)+M_{0}}{D_{0}}}\qquad {\mbox{ Calc the numeric value }}$

Step 4. Substitute $X_{0}$ with $N_{1}+{\frac {1}{X_{1}}}$

$X_{0LowerBound}\quad <\quad X_{0}\quad <\quad X_{0UpperBound}$

after substitution

$X_{0LowerBound}\quad <\quad N_{1}+{\frac {1}{X_{1}}}\quad <\quad X_{0UpperBound}$

Since ${\frac {1}{X_{1}}}$ is less than one, we can determine $N_{1}$ from the numeric values of $X_{0LowerBound}$ and $X_{0UpperBound}$ because $N_{1}\ \in \ \mathbb {Z^{+}}$ (ie. $N_{1}$ is an integer). Determine the numeric value of $N_{1}$.

Step 5. Once we know the value of $N_{1}$, we can rework the equation for $X_{1}$

$X_{0}=N_{1}+{\frac {1}{X_{1}}}\quad \longrightarrow \quad X_{1}={\frac {1}{X_{0}-N_{1}}}$

$X_{1}={\frac {1}{{\frac {{\sqrt {r}}+M_{0}}{D_{0}}}-N_{1}}}={\frac {1}{\frac {{\sqrt {r}}+M_{0}-D_{0}\,N_{1}}{D_{0}}}}={\frac {D_{0}}{{\sqrt {r}}+M_{0}-D_{0}\,N_{1}}}$

$X_{1}={\frac {D_{0}}{{\sqrt {r}}+(M_{0}-D_{0}\,N_{1})}}\times {\frac {{\sqrt {r}}-(M_{0}-D_{0}\,N_{1})}{{\sqrt {r}}-(M_{0}-D_{0}\,N_{1})}}={\frac {D_{0}({\sqrt {r}}-(M_{0}-D_{0}\,N_{1}))}{r-(M_{0}-D_{0}\,N_{1})^{2}}}$

$X_{1}={\frac {{\sqrt {r}}-(M_{0}-D_{0}\,N_{1})}{\frac {r-(M_{0}-D_{0}\,N_{1})^{2}}{D_{0}}}}={\frac {{\sqrt {r}}+M_{1}}{D_{1}}}$

${\mbox{ where }}\quad M_{1}=-(M_{0}-D_{0}\,N_{1})$

and

${\mbox{ where }}\quad D_{1}={\frac {r-(M_{0}-D_{0}\,N_{1})^{2}}{D_{0}}}$

Step 6. Write $X_{1}$ in terms of ${\sqrt {r}}$.

$X_{1}={\frac {{\sqrt {r}}+M_{1}}{D_{1}}}$

$\downarrow$

${\mbox{ substitute }}\quad {\sqrt {r}}\quad {\mbox{ with }}\quad L\quad {\mbox{ or }}\quad U$

$\downarrow$

$X_{1LowerBound}={\frac {{\sqrt {r}}+M_{1}}{D_{1}}}={\frac {L+M_{1}}{D_{1}}}={\frac {N_{0}+M_{1}}{D_{1}}}\qquad {\mbox{ Calc the numeric value }}$

and

$X_{1UpperBound}={\frac {{\sqrt {r}}+M_{1}}{D_{1}}}={\frac {U+M_{1}}{D_{1}}}={\frac {(N_{0}+1)+M_{1}}{D_{1}}}\qquad {\mbox{ Calc the numeric value }}$

Step 7. Substitute $X_{1}$ with $N_{2}+{\frac {1}{X_{2}}}$

$X_{1LowerBound}\quad <\quad X_{1}\quad <\quad X_{1UpperBound}$

after substitution

$X_{1LowerBound}\quad <\quad N_{2}+{\frac {1}{X_{2}}}\quad <\quad X_{1UpperBound}$

Since ${\frac {1}{X_{2}}}$ is less than one, we can determine $N_{2}$ from the numeric values of $X_{1LowerBound}$ and $X_{1UpperBound}$ because $N_{2}\ \in \ \mathbb {Z^{+}}$ (ie. $N_{2}$ is an integer). Determine the numeric value of $N_{2}$.

Step 8. Once we know the value of $N_{2}$, we can rework the equation for $X_{2}$

$X_{1}=N_{2}+{\frac {1}{X_{2}}}\quad \longrightarrow \quad X_{2}={\frac {1}{X_{1}-N_{2}}}$

$X_{2}={\frac {1}{{\frac {{\sqrt {r}}+M_{1}}{D_{1}}}-N_{2}}}={\frac {1}{\frac {{\sqrt {r}}+M_{1}-D_{1}\,N_{2}}{D_{1}}}}={\frac {D_{1}}{{\sqrt {r}}+M_{1}-D_{1}\,N_{2}}}$

$X_{2}={\frac {D_{1}}{{\sqrt {r}}+(M_{1}-D_{1}\,N_{2})}}\times {\frac {{\sqrt {r}}-(M_{1}-D_{1}\,N_{2})}{{\sqrt {r}}-(M_{1}-D_{1}\,N_{2})}}={\frac {D_{1}({\sqrt {r}}-(M_{1}-D_{1}\,N_{2}))}{r-(M_{1}-D_{1}\,N_{2})^{2}}}$

$X_{2}={\frac {{\sqrt {r}}-(M_{1}-D_{1}\,N_{2})}{\frac {r-(M_{1}-D_{1}\,N_{2})^{2}}{D_{1}}}}={\frac {{\sqrt {r}}+M_{2}}{D_{2}}}$

${\mbox{ where }}\quad M_{2}=-(M_{1}-D_{1}\,N_{2})$

and

${\mbox{ where }}\quad D_{2}={\frac {r-(M_{1}-D_{1}\,N_{2})^{2}}{D_{1}}}$

Step 9. Repeat step 6 , 7 and 8 .

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.