Methods of computing square roots

Typically the number ${\displaystyle S}$ is expressed in scientific notation as ${\displaystyle a\times 10^{2n}}$ where ${\displaystyle 1\leq a<100}$ and n is an integer, and the range of possible square roots is ${\displaystyle {\sqrt {a}}\times 10^{n}}$ where ${\displaystyle 1\leq {\sqrt {a}}<10}$.

Scalar methods divide the range into intervals, and the estimate in each interval is represented by a single scalar number. If the range is considered as a single interval, the arithmetic mean (5.5) or geometric mean (${\displaystyle {\sqrt {10}}\approx 3.16}$) times ${\displaystyle 10^{n}}$ are plausible estimates. The absolute and relative error for these will differ. In general, a single scalar will be very inaccurate. Better estimates divide the range into two or more intervals, but scalar estimates have inherently low accuracy.

For two intervals, divided geometrically, the square root ${\displaystyle {\sqrt {S}}={\sqrt {a}}\times 10^{n}}$ can be estimated as[Note 2]

${\displaystyle {\sqrt {S}}\approx {\begin{cases}2\cdot 10^{n}&{\text{if }}a<10,\\6\cdot 10^{n}&{\text{if }}a\geq 10.\end{cases}}}$

This estimate has maximum absolute error of ${\displaystyle 4\cdot 10^{n}}$ at a = 100, and maximum relative error of 100% at a = 1.

For example, for ${\displaystyle S=125348}$ factored as ${\displaystyle 12.5348\times 10^{4}}$, the estimate is ${\displaystyle {\sqrt {S}}\approx 6\cdot 10^{2}=600}$. ${\displaystyle {\sqrt {125348}}=354.0}$, an absolute error of 246 and relative error of almost 70%.

A better estimate, and the standard method used, is a linear approximation to the function ${\displaystyle y=x^{2}}$ over a small arc. If, as above, powers of the base are factored out of the number ${\displaystyle S}$ and the interval reduced to ${\displaystyle [1,100]}$, a secant line spanning the arc, or a tangent line somewhere along the arc may be used as the approximation, but a least-squares regression line intersecting the arc will be more accurate.

A least-squares regression line minimizes the average difference between the estimate and the value of the function. Its equation is ${\displaystyle y=8.7x-10}$. Reordering, ${\displaystyle x=0.115y+1.15}$. Rounding the coefficients for ease of computation,

${\displaystyle {\sqrt {S}}\approx (a/10+1.2)\cdot 10^{n}}$

That is the best estimate on average that can be achieved with a single piece linear approximation of the function y=x² in the interval ${\displaystyle [1,100]}$. It has a maximum absolute error of 1.2 at a=100, and maximum relative error of 30% at S=1 and 10.[Note 3]

To divide by 10, subtract one from the exponent of ${\displaystyle a}$, or figuratively move the decimal point one digit to the left. For this formulation, any additive constant 1 plus a small increment will make a satisfactory estimate so remembering the exact number isn't a burden. The approximation (rounded or not) using a single line spanning the range ${\displaystyle [1,100]}$ is less than one significant digit of precision; the relative error is greater than 1/2², so less than 2 bits of information are provided. The accuracy is severely limited because the range is two orders of magnitude, quite large for this kind of estimation.

A much better estimate can be obtained by a piece-wise linear approximation: multiple line segments, each approximating some subarc of the original. The more line segments used, the better the approximation. The most common way is to use tangent lines; the critical choices are how to divide the arc and where to place the tangent points. An efficacious way to divide the arc from y=1 to y=100 is geometrically: for two intervals, the bounds of the intervals are the square root of the bounds of the original interval, 1*100, i.e. [1,²√100] and [²√100,100]. For three intervals, the bounds are the cube roots of 100: [1,³√100], [³√100,(³√100)²], and [(³√100)²,100], etc. For two intervals, ²√100 = 10, a very convenient number. Tangent lines are easy to derive, and are located at x = √1*√10 and x = √10*√10. Their equations are: ${\displaystyle y=3.56x-3.16}$ and ${\displaystyle y=11.2x-31.6}$. Inverting, the square roots are: ${\displaystyle x=0.28y+0.89}$ and ${\displaystyle x=.089y+2.8}$. Thus for ${\displaystyle S=a\cdot 10^{2n}}$:

${\displaystyle {\sqrt {S}}\approx {\begin{cases}(0.28a+0.89)\cdot 10^{n}&{\text{if }}a<10,\\(.089a+2.8)\cdot 10^{n}&{\text{if }}a\geq 10.\end{cases}}}$

The maximum absolute errors occur at the high points of the intervals, at a=10 and 100, and are 0.54 and 1.7 respectively. The maximum relative errors are at the endpoints of the intervals, at a=1, 10 and 100, and are 17% in both cases. 17% or 0.17 is larger than 1/10, so the method yields less than a decimal digit of accuracy.

In some cases, hyperbolic estimates may be efficacious, because a hyperbola is also a convex curve and may lie along an arc of Y = x² better than a line. Hyperbolic estimates are more computationally complex, because they necessarily require a floating division. A near-optimal hyperbolic approximation to x² on the interval ${\displaystyle [1,100]}$ is y=190/(10-x)-20. Transposing, the square root is x = -190/(y+20)+10. Thus for ${\displaystyle S=a\cdot 10^{2n}}$:

${\displaystyle {\sqrt {S}}\approx \left({\frac {-190}{a+20}}+10\right)\cdot 10^{n}}$

The floating division need be accurate to only one decimal digit, because the estimate overall is only that accurate, and can be done mentally. A hyperbolic estimate is better on average than scalar or linear estimates. It has maximum absolute error of 1.58 at 100 and maximum relative error of 16.0% at 10. For the worst case at a=10, the estimate is 3.67. If one starts with 10 and applies Newton-Raphson iterations straight away, two iterations will be required, yielding 3.66, before the accuracy of the hyperbolic estimate is exceeded. For a more typical case like 75, the hyperbolic estimate is 8.00, and 5 Newton-Raphson iterations starting at 75 would be required to obtain a more accurate result.

A method analogous to piece-wise linear approximation but using only arithmetic instead of algebraic equations, uses the multiplication tables in reverse: the square root of a number between 1 and 100 is between 1 and 10, so if we know 25 is a perfect square (5 × 5), and 36 is a perfect square (6 × 6), then the square root of a number greater than or equal to 25 but less than 36, begins with a 5. Similarly for numbers between other squares. This method will yield a correct first digit, but it is not accurate to one digit: the first digit of the square root of 35 for example, is 5, but the square root of 35 is almost 6.

A better way is to the divide the range into intervals half way between the squares. So any number between 25 and half way to 36, which is 30.5, estimate 5; any number greater than 30.5 up to 36, estimate 6.[Note 4] The procedure only requires a little arithmetic to find a boundary number in the middle of two products from the multiplication table. Here is a reference table of those boundaries:

a	nearest square	${\displaystyle k={\sqrt {a}}}$ est.
1 to 2.5	1 (= 12)	1
2.5 to 6.5	4 (= 22)	2
6.5 to 12.5	9 (= 32)	3
12.5 to 20.5	16 (= 42)	4
20.5 to 30.5	25 (= 52)	5
30.5 to 42.5	36 (= 62)	6
42.5 to 56.5	49 (= 72)	7
56.5 to 72.5	64 (= 82)	8
72.5 to 90.5	81 (= 92)	9
90.5 to 100	100 (= 102)	10

The final operation is to multiply the estimate k by the power of ten divided by 2, so for ${\displaystyle S=a\cdot 10^{2n}}$,

${\displaystyle {\sqrt {S}}\approx k\cdot 10^{n}}$

The method implicitly yields one significant digit of accuracy, since it rounds to the best first digit.

The method can be extended 3 significant digits in most cases, by interpolating between the nearest squares bounding the operand. If ${\displaystyle k^{2}\leq a<(k+1)^{2}}$, then ${\displaystyle {\sqrt {a}}}$ is approximately k plus a fraction, the difference between a and k² divided by the difference between the two squares:

${\displaystyle {\sqrt {a}}\approx k+R}$ where ${\displaystyle R={\frac {(a-k^{2})}{(k+1)^{2}-k^{2}}}}$

The final operation, as above, is to multiply the result by the power of ten divided by 2;

${\displaystyle {\sqrt {S}}={\sqrt {a}}\cdot 10^{n}\approx (k+R)\cdot 10^{n}}$

k is a decimal digit and R is a fraction that must be converted to decimal. It usually has only a single digit in the numerator, and one or two digits in the denominator, so the conversion to decimal can be done mentally.

Example: find the square root of 75. 75 = 75 × 10^{2 · 0}, so a is 75 and n is 0. From the multiplication tables, the square root of the mantissa must be 8 point something because 8 × 8 is 64, but 9 × 9 is 81, too big, so k is 8; something is the decimal representation of R. The fraction R is 75 - k² = 11, the numerator, and 81 - k² = 17, the denominator. 11/17 is a little less than 12/18, which is 2/3s or .67, so guess .66 (it's ok to guess here, the error is very small). So the estimate is 8 + .66 = 8.66. √75 to three significant digits is 8.66, so the estimate is good to 3 significant digits. Not all such estimates using this method will be so accurate, but they will be close.

When working in the binary numeral system (as computers do internally), by expressing ${\displaystyle S}$ as ${\displaystyle a\times 2^{2n}}$ where ${\displaystyle 0.1_{2}\leq a<10_{2}}$, the square root ${\displaystyle {\sqrt {S}}={\sqrt {a}}\times 2^{n}}$ can be estimated as

${\displaystyle {\sqrt {S}}\approx (0.485+0.485\cdot a)\cdot 2^{n}}$

which is the least-squares regression line to 3 significant digit coefficients. ${\displaystyle {\sqrt {a}}}$ has maximum absolute error of 0.0408 at ${\displaystyle a}$=2, and maximum relative error of 3.0% at ${\displaystyle a}$=1. A computationally convenient rounded estimate (because the coefficients are powers of 2) is:

${\displaystyle {\sqrt {S}}\approx (0.5+0.5\cdot a)\cdot 2^{n}}$[Note 5]

which has maximum absolute error of 0.086 at 2 and maximum relative error of 6.1% at ${\displaystyle a}$=0.5 and ${\displaystyle a}$=2.0.

For ${\displaystyle S=125348=1\;1110\;1001\;1010\;0100_{2}=1.1110\;1001\;1010\;0100_{2}\times 2^{16}\,}$, the binary approximation gives ${\displaystyle {\sqrt {S}}\approx (0.5+0.5\cdot a)\cdot 2^{8}=1.0111\;0100\;1101\;0010_{2}\cdot 1\;0000\;0000_{2}=1.456\cdot 256=372.8}$. ${\displaystyle {\sqrt {125348}}=354.0}$, so the estimate has an absolute error of 19 and relative error of 5.3%. The relative error is a little less than 1/2⁴, so the estimate is good to 4+ bits.

An estimate for ${\displaystyle a}$ good to 8 bits can be obtained by table lookup on the high 8 bits of ${\displaystyle a}$, remembering that the high bit is implicit in most floating point representations, and the bottom bit of the 8 should be rounded. The table is 256 bytes of precomputed 8-bit square root values. For example, for the index 11101101₂ representing 1.8515625₁₀, the entry is 10101110₂ representing 1.359375₁₀, the square root of 1.8515625₁₀ to 8 bit precision (2+ decimal digits).

To calculate √S, where S = 125348, to six significant figures, use the rough estimation method above to get

${\displaystyle {\begin{aligned}{\begin{array}{rlll}x_{0}&=6\cdot 10^{2}&&=600.000\\[0.3em]x_{1}&={\frac {1}{2}}\left(x_{0}+{\frac {S}{x_{0}}}\right)&={\frac {1}{2}}\left(600.000+{\frac {125348}{600.000}}\right)&=404.457\\[0.3em]x_{2}&={\frac {1}{2}}\left(x_{1}+{\frac {S}{x_{1}}}\right)&={\frac {1}{2}}\left(404.457+{\frac {125348}{404.457}}\right)&=357.187\\[0.3em]x_{3}&={\frac {1}{2}}\left(x_{2}+{\frac {S}{x_{2}}}\right)&={\frac {1}{2}}\left(357.187+{\frac {125348}{357.187}}\right)&=354.059\\[0.3em]x_{4}&={\frac {1}{2}}\left(x_{3}+{\frac {S}{x_{3}}}\right)&={\frac {1}{2}}\left(354.059+{\frac {125348}{354.059}}\right)&=354.045\\[0.3em]x_{5}&={\frac {1}{2}}\left(x_{4}+{\frac {S}{x_{4}}}\right)&={\frac {1}{2}}\left(354.045+{\frac {125348}{354.045}}\right)&=354.045\end{array}}\end{aligned}}}$

Therefore, √125348 ≈ 354.045.

Semilog graphs comparing the speed of convergence of Heron's method to find the square root of 100 for different initial guesses. Negative guesses converge to the negative root, positive guesses to the positive root. Note that values closer to the root converge faster, and all approximations are overestimates. In the SVG file, hover over a graph to display its points.

Suppose that x₀ > 0 and S > 0. Then for any natural number n, x_n > 0. Let the relative error in x_n be defined by

${\displaystyle \varepsilon _{n}={\frac {x_{n}}{\sqrt {S}}}-1>-1}$

and thus

${\displaystyle x_{n}={\sqrt {S}}\cdot (1+\varepsilon _{n}).}$

Then it can be shown that

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

And thus that

${\displaystyle \varepsilon _{n+2}\leq \min \left\{{\frac {\varepsilon _{n+1}^{2}}{2}},{\frac {\varepsilon _{n+1}}{2}}\right\}}$

and consequently that convergence is assured, and quadratic.

If using the rough estimate above with the Babylonian method, then the least accurate cases in ascending order are as follows:

${\displaystyle {\begin{aligned}S&=1;&x_{0}&=2;&x_{1}&=1.250;&\varepsilon _{1}&=0.250.\\S&=10;&x_{0}&=2;&x_{1}&=3.500;&\varepsilon _{1}&<0.107.\\S&=10;&x_{0}&=6;&x_{1}&=3.833;&\varepsilon _{1}&<0.213.\\S&=100;&x_{0}&=6;&x_{1}&=11.333;&\varepsilon _{1}&<0.134.\end{aligned}}}$

Thus in any case,

${\displaystyle \varepsilon _{1}\leq 2^{-2}.\,}$

${\displaystyle \varepsilon _{2}<2^{-5}<10^{-1}.\,}$

${\displaystyle \varepsilon _{3}<2^{-11}<10^{-3}.\,}$

${\displaystyle \varepsilon _{4}<2^{-23}<10^{-6}.\,}$

${\displaystyle \varepsilon _{5}<2^{-47}<10^{-14}.\,}$

${\displaystyle \varepsilon _{6}<2^{-95}<10^{-28}.\,}$

${\displaystyle \varepsilon _{7}<2^{-191}<10^{-57}.\,}$

${\displaystyle \varepsilon _{8}<2^{-383}<10^{-115}.\,}$

Rounding errors will slow the convergence. It is recommended to keep at least one extra digit beyond the desired accuracy of the x_n being calculated to minimize round off error.

Using the same example as given with the Babylonian method, let ${\displaystyle S=125348.}$ Then, the first iteration gives

${\displaystyle {\begin{aligned}x_{0}&=600\\a_{1}&={\frac {125348-600^{2}}{2\times 600}}&&=&-195.543\\b_{1}&=600+(-195.543)&&=&404.456\\x_{1}&=404.456-{\frac {(-195.543)^{2}}{2\times 404.456}}&&=&357.186\end{aligned}}}$

Likewise the second iteration gives

${\displaystyle {\begin{aligned}a_{2}&={\frac {125348-357.186^{2}}{2\times 357.186}}&&=&-3.126\\b_{2}&=357.186+(-3.126)&&=&354.060\\x_{2}&=354.06-{\frac {(-3.1269)^{2}}{2\times 354.06}}&&=&354.046\end{aligned}}}$

First, consider the case of finding the square root of a number Z, that is the square of a two-digit number XY, where X is the tens digit and Y is the units digit. Specifically:

${\displaystyle Z=(10X+Y)^{2}=100X^{2}+20XY+Y^{2}}$

Now using the digit-by-digit algorithm, we first determine the value of X. X is the largest digit such that X² is less than or equal to Z from which we removed the two rightmost digits.

In the next iteration, we pair the digits, multiply X by 2, and place it in the tenth's place while we try to figure out what the value of Y is.

Since this is a simple case where the answer is a perfect square root XY, the algorithm stops here.

The same idea can be extended to any arbitrary square root computation next. Suppose we are able to find the square root of N by expressing it as a sum of n positive numbers such that

${\displaystyle N=(a_{1}+a_{2}+a_{3}+\dotsb +a_{n})^{2}}$.

By repeatedly applying the basic identity

${\displaystyle (x+y)^{2}=x^{2}+2xy+y^{2},}$

the right-hand-side term can be expanded as

${\displaystyle {\begin{aligned}&(a_{1}+a_{2}+a_{3}+\dotsb +a_{n})^{2}\\=&\,a_{1}^{2}+2a_{1}a_{2}+a_{2}^{2}+2(a_{1}+a_{2})a_{3}+a_{3}^{2}+\dotsb +a_{n-1}^{2}+2\left(\sum _{i=1}^{n-1}a_{i}\right)a_{n}+a_{n}^{2}\\=&\,a_{1}^{2}+[2a_{1}+a_{2}]a_{2}+[2(a_{1}+a_{2})+a_{3}]a_{3}+\dotsb +\left[2\left(\sum _{i=1}^{n-1}a_{i}\right)+a_{n}\right]a_{n}.\end{aligned}}}$

This expression allows us to find the square root by sequentially guessing the values of ${\displaystyle a_{i}}$s. Suppose that the numbers ${\displaystyle a_{1},\ldots ,a_{m-1}}$ have already been guessed, then the m-th term of the right-hand-side of above summation is given by ${\displaystyle Y_{m}=[2P_{m-1}+a_{m}]a_{m},}$ where ${\displaystyle P_{m-1}=\sum _{i=1}^{m-1}a_{i}}$ is the approximate square root found so far. Now each new guess ${\displaystyle a_{m}}$ should satisfy the recursion

${\displaystyle X_{m}=X_{m-1}-Y_{m},}$

such that ${\displaystyle X_{m}\geq 0}$ for all ${\displaystyle 1\leq m\leq n,}$ with initialization ${\displaystyle X_{0}=N.}$ When ${\displaystyle X_{n}=0,}$ the exact square root has been found; if not, then the sum of ${\displaystyle a_{i}}$s gives a suitable approximation of the square root, with ${\displaystyle X_{n}}$ being the approximation error.

For example, in the decimal number system we have

${\displaystyle N=(a_{1}\cdot 10^{n-1}+a_{2}\cdot 10^{n-2}+\cdots +a_{n-1}\cdot 10+a_{n})^{2},}$

where ${\displaystyle 10^{n-i}}$ are place holders and the coefficients ${\displaystyle a_{i}\in \{0,1,2,\ldots ,9\}}$. At any m-th stage of the square root calculation, the approximate root found so far, ${\displaystyle P_{m-1}}$ and the summation term ${\displaystyle Y_{m}}$ are given by

${\displaystyle P_{m-1}=\sum _{i=1}^{m-1}a_{i}\cdot 10^{n-i}=10^{n-m+1}\sum _{i=1}^{m-1}a_{i}\cdot 10^{m-i-1},}$

${\displaystyle Y_{m}=[2P_{m-1}+a_{m}\cdot 10^{n-m}]a_{m}\cdot 10^{n-m}=\left[20\sum _{i=1}^{m-1}a_{i}\cdot 10^{m-i-1}+a_{m}\right]a_{m}\cdot 10^{2(n-m)}.}$

Here since the place value of ${\displaystyle Y_{m}}$ is an even power of 10, we only need to work with the pair of most significant digits of the remaining term ${\displaystyle X_{m-1}}$ at any m-th stage. The section below codifies this procedure.

It is obvious that a similar method can be used to compute the square root in number systems other than the decimal number system. For instance, finding the digit-by-digit square root in the binary number system is quite efficient since the value of ${\displaystyle a_{i}}$ is searched from a smaller set of binary digits {0,1}. This makes the computation faster since at each stage the value of ${\displaystyle Y_{m}}$ is either ${\displaystyle Y_{m}=0}$ for ${\displaystyle a_{m}=0}$ or ${\displaystyle Y_{m}=2P_{m-1}+1}$ for ${\displaystyle a_{m}=1}$. The fact that we have only two possible options for ${\displaystyle a_{m}}$ also makes the process of deciding the value of ${\displaystyle a_{m}}$ at m-th stage of calculation easier. This is because we only need to check if ${\displaystyle Y_{m}\leq X_{m-1}}$ for ${\displaystyle a_{m}=1.}$ If this condition is satisfied, then we take ${\displaystyle a_{m}=1}$; if not then ${\displaystyle a_{m}=0.}$ Also, the fact that multiplication by 2 is done by left bit-shifts helps in the computation.

Write the original number in decimal form. The numbers are written similar to the long division algorithm, and, as in long division, the root will be written on the line above. Now separate the digits into pairs, starting from the decimal point and going both left and right. The decimal point of the root will be above the decimal point of the square. One digit of the root will appear above each pair of digits of the square.

Beginning with the left-most pair of digits, do the following procedure for each pair:

1. Starting on the left, bring down the most significant (leftmost) pair of digits not yet used (if all the digits have been used, write "00") and write them to the right of the remainder from the previous step (on the first step, there will be no remainder). In other words, multiply the remainder by 100 and add the two digits. This will be the current value c.
2. Find p, y and x, as follows:
   • Let p be the part of the root found so far, ignoring any decimal point. (For the first step, p = 0.)
   • Determine the greatest digit x such that ${\displaystyle x(20p+x)\leq c}$. We will use a new variable y = x(20p + x).
     • Note: 20p + x is simply twice p, with the digit x appended to the right.
     • Note: x can be found by guessing what c/(20·p) is and doing a trial calculation of y, then adjusting x upward or downward as necessary.
   • Place the digit ${\displaystyle x}$ as the next digit of the root, i.e., above the two digits of the square you just brought down. Thus the next p will be the old p times 10 plus x.
3. Subtract y from c to form a new remainder.
4. If the remainder is zero and there are no more digits to bring down, then the algorithm has terminated. Otherwise go back to step 1 for another iteration.

Find the square root of 152.2756.

          1  2. 3  4 
       /
     \/  01 52.27 56

         01                   1*1 <= 1 < 2*2                 x=1
         01                     y = x*x = 1*1 = 1
         00 52                22*2 <= 52 < 23*3              x=2
         00 44                  y = (20+x)*x = 22*2 = 44
            08 27             243*3 <= 827 < 244*4           x=3
            07 29               y = (240+x)*x = 243*3 = 729
               98 56          2464*4 <= 9856 < 2465*5        x=4
               98 56            y = (2460+x)*x = 2464*4 = 9856
               00 00          Algorithm terminates: Answer=12.34

This section uses the formalism from the digit-by-digit calculation section above, with the slight variation that we let ${\displaystyle N^{2}=(a_{n}+\dotsb +a_{0})^{2}}$, with each ${\displaystyle a_{m}=2^{m}}$ or ${\displaystyle a_{m}=0}$.
We iterate all ${\displaystyle 2^{m}}$, from ${\displaystyle 2^{n}}$ down to ${\displaystyle 2^{0}}$, and build up an approximate solution ${\displaystyle P_{m}=a_{n}+a_{n-1}+\ldots +a_{m}}$, the sum of all ${\displaystyle a_{i}}$ for which we have determined the value.
To determine if ${\displaystyle a_{m}}$ equals ${\displaystyle 2^{m}}$ or ${\displaystyle 0}$, we let ${\displaystyle P_{m}=P_{m+1}+2^{m}}$. If ${\displaystyle P_{m}^{2}\leq N^{2}}$ (i.e. the square of our approximate solution including ${\displaystyle 2^{m}}$ does not exceed the target square) then ${\displaystyle a_{m}=2^{m}}$, otherwise ${\displaystyle a_{m}=0}$ and ${\displaystyle P_{m}=P_{m+1}}$.
To avoid squaring ${\displaystyle P_{m}}$ in each step, we store the difference ${\displaystyle X_{m}=N^{2}-P_{m}^{2}}$ and incrementally update it by setting ${\displaystyle X_{m}=X_{m+1}-Y_{m}}$ with ${\displaystyle Y_{m}=P_{m}^{2}-P_{m+1}^{2}=2P_{m+1}a_{m}+a_{m}^{2}}$.
Initially, we set ${\displaystyle a_{n}=P_{n}=2^{n}}$ for the largest ${\displaystyle n}$ with ${\displaystyle (2^{n})^{2}=4^{n}\leq N^{2}}$.

As an extra optimization, we store ${\displaystyle P_{m+1}2^{m+1}}$ and ${\displaystyle (2^{m})^{2}}$, the two terms of ${\displaystyle Y_{m}}$ in case that ${\displaystyle a_{m}}$ is nonzero, in separate variables ${\displaystyle c_{m}}$, ${\displaystyle d_{m}}$:

${\displaystyle c_{m}=P_{m+1}2^{m+1}}$

${\displaystyle d_{m}=(2^{m})^{2}}$

${\displaystyle Y_{m}={\begin{cases}c_{m}+d_{m}&{\text{if }}a_{m}=2^{m}\\0&{\text{if }}a_{m}=0\end{cases}}}$

${\displaystyle c_{m}}$ and ${\displaystyle d_{m}}$ can be efficiently updated in each step:

${\displaystyle c_{m-1}=P_{m}2^{m}=(P_{m+1}+a_{m})2^{m}=P_{m+1}2^{m}+a_{m}2^{m}={\begin{cases}c_{m}/2+d_{m}&{\text{if }}a_{m}=2^{m}\\c_{m}/2&{\text{if }}a_{m}=0\end{cases}}}$

${\displaystyle d_{m-1}={\frac {d_{m}}{4}}}$

Note that:

${\displaystyle c_{-1}=P_{0}2^{0}=P_{0}=N}$, which is the final result returned in the function below.

An implementation of this algorithm in C:[5]

int32_t isqrt(int32_t n) {
    assert(("sqrt input should be non-negative", n > 0));

    // Xₙ₊₁
    int32_t x = n;

    // cₙ
    int32_t c = 0;

    // dₙ which starts at the highest power of four <= n
    int32_t d = 1 << 30; // The second-to-top bit is set.
                         // Same as ((unsigned) INT32_MAX + 1) / 2.
    while (d > n)
        d >>= 2;

    // for dₙ … d₀
    while (d != 0) {
        if (x >= c + d) {      // if Xₘ₊₁ ≥ Yₘ then aₘ = 2ᵐ
            x -= c + d;        // Xₘ = Xₘ₊₁ - Yₘ
            c = (c >> 1) + d;  // cₘ₋₁ = cₘ/2 + dₘ (aₘ is 2ᵐ)
        }
        else {
            c >>= 1;           // cₘ₋₁ = cₘ/2      (aₘ is 0)
        }
        d >>= 2;               // dₘ₋₁ = dₘ/4
    }
    return c;                  // c₋₁
}

Faster algorithms, in binary and decimal or any other base, can be realized by using lookup tables—in effect trading more storage space for reduced run time.[6]

Some computers use Goldschmidt's algorithm to simultaneously calculate ${\displaystyle {\sqrt {S}}}$ and ${\displaystyle 1/{\sqrt {S}}}$. Goldschmidt's algorithm finds ${\displaystyle {\sqrt {S}}}$ faster than Newton-Raphson iteration on a computer with a fused multiply–add instruction and either a pipelined floating point unit or two independent floating-point units.[10]

The first way of writing Goldschmidt's algorithm begins

${\displaystyle b_{0}=S}$

${\displaystyle Y_{0}\approx 1/{\sqrt {S}}}$ (typically using a table lookup)

${\displaystyle y_{0}=Y_{0}}$

${\displaystyle x_{0}=Sy_{0}}$

and iterates

${\displaystyle b_{n+1}=b_{n}Y_{n}^{2}}$

${\displaystyle Y_{n+1}=(3-b_{n+1})/2}$

${\displaystyle x_{n+1}=x_{n}Y_{n+1}}$

${\displaystyle y_{n+1}=y_{n}Y_{n+1}}$

until ${\displaystyle b_{i}}$ is sufficiently close to 1, or a fixed number of iterations. The iterations converge to

${\displaystyle \lim _{n\to \infty }x_{n}={\sqrt {S}}}$, and

${\displaystyle \lim _{n\to \infty }y_{n}=1/{\sqrt {S}}}$.

Note that it is possible to omit either ${\displaystyle x_{n}}$ and ${\displaystyle y_{n}}$ from the computation, and if both are desired then ${\displaystyle x_{n}=Sy_{n}}$ may be used at the end rather than computing it through in each iteration.

A second form, using fused multiply-add operations, begins

${\displaystyle y_{0}\approx 1/{\sqrt {S}}}$ (typically using a table lookup)

${\displaystyle x_{0}=Sy_{0}}$

${\displaystyle h_{0}=y_{0}/2}$

and iterates

${\displaystyle r_{n}=0.5-x_{n}h_{n}}$

${\displaystyle x_{n+1}=x_{n}+x_{n}r_{n}}$

${\displaystyle h_{n+1}=h_{n}+h_{n}r_{n}}$

until ${\displaystyle r_{i}}$ is sufficiently close to 0, or a fixed number of iterations. This converges to

${\displaystyle \lim _{n\to \infty }x_{n}={\sqrt {S}}}$, and

${\displaystyle \lim _{n\to \infty }2h_{n}=1/{\sqrt {S}}}$.