View source for Methods of computing square roots

{{short description|Algorithms for calculating square roots}}
{{multiple issues|
{{Original research|date=January 2012}}
{{technical|date=September 2012}}
}}

'''Methods of computing square roots''' are [[algorithm]]s for approximating the non-negative [[square root]] <math>\sqrt{S}</math> of a positive [[real number]] <math>S</math>.
Since all square roots of [[natural number]]s, other than of [[square number|perfect square]]s, are [[Irrational number|irrational]],{{sfn|Jackson|2011}}
square roots can usually only be computed to some finite precision: these methods typically construct a series of increasingly accurate [[Numerical approximation|approximations]].

Most square root computation methods are iterative: after choosing a suitable [[#Initial estimate|initial estimate]] of <math>\sqrt{S}</math>, an iterative refinement is performed until some termination criterion is met.
One refinement scheme is [[#Heron's method|Heron's method]], a special case of [[Newton's method]].
If division is much more costly than multiplication, it may be preferable to compute the [[#Iterative methods for reciprocal square roots|inverse square root]] instead.

Other methods are available to compute the square root [[#Digit-by-digit calculation|digit by digit]], or using [[#Taylor series|Taylor series]].
Rational approximations of square roots may be calculated using [[#Continued fraction expansion|continued fraction expansions]].

The method employed depends on the needed accuracy, and the available tools and computational power. The methods may be roughly classified as those suitable for mental calculation, those usually requiring at least paper and pencil, and those which are implemented as programs to be executed on a digital electronic computer or other computing device. Algorithms may take into account convergence (how many iterations are required to achieve a specified precision), computational complexity of individual operations (i.e. division) or iterations, and error propagation (the accuracy of the final result).

A few methods like paper-and-pencil synthetic division and series expansion, do not require a starting value. In some applications, an [[integer square root]] is required, which is the square root rounded or truncated to the nearest integer (a modified procedure may be employed in this case).

==History==
Procedures for finding square roots (particularly the [[square root of 2]]) have been known since at least the period of ancient Babylon in the 17th century BCE. 
[[Babylonian mathematics|Babylonian mathematicians]] calculated the square root of 2 to three sexagesimal "digits" after the 1, but it is not known exactly how. They knew how to approximate a hypotenuse using
<math display="block">\sqrt{a^2+b^2}\approx a+\frac{b^2}{2a}</math>
(giving for example <math>\frac{41}{60}+\frac{15}{3600}</math> for the diagonal of a gate whose height is <math>\frac{40}{60}</math> rods and whose width is <math>\frac{10}{60}</math> rods) and they may have used a similar approach for finding the approximation of <math>\sqrt 2.</math>{{sfn|Fowler|Robson|1998}}

[[Heron's method]] from first century Egypt was the first ascertainable algorithm for computing square root.{{sfn|Heath|1921}}

Modern analytic methods began to be developed after introduction of the [[Arabic numeral]] system to western Europe in the early Renaissance.{{cn|date=December 2023}}

Today, nearly all computing devices have a fast and accurate square root function, either as a programming language construct, a compiler intrinsic or library function, or as a hardware operator, based on one of the described procedures.

==Initial estimate==
Many iterative square root algorithms require an initial [[seed value]]. The seed must be a non-zero positive number; it should be between 1 and <math>S</math>, the number whose square root is desired, because the square root must be in that range. If the seed is far away from the root, the algorithm will require more iterations. If one initializes with <math>x_0=1</math> (or <math>S</math>), then approximately <math> \tfrac12 \vert \log_2 S \vert </math> iterations will be wasted just getting the order of magnitude of the root. It is therefore useful to have a rough estimate, which may have limited accuracy but is easy to calculate. In general, the better the initial estimate, the faster the convergence.  For Newton's method, a seed somewhat larger than the root will converge slightly faster than a seed somewhat smaller than the root.

In general, an estimate is pursuant to an arbitrary interval known to contain the root (such as <math>[x_0,S/x_0]</math>). The estimate is a specific value of a functional approximation to <math>f(x)=\sqrt{x}</math> over the interval. Obtaining a better estimate involves either obtaining tighter bounds on the interval, or finding a better functional approximation to <math>f(x)</math>. The latter usually means using a higher order polynomial in the approximation, though not all approximations are polynomial. Common methods of estimating include scalar, linear, hyperbolic and logarithmic. A decimal base is usually used for mental or paper-and-pencil estimating. A binary base is more suitable for computer estimates.  In estimating, the exponent and [[mantissa (logarithm)|mantissa]] are usually treated separately, as the number would be expressed in scientific notation.

===Decimal estimates===
Typically the number <math>S</math> is expressed in [[scientific notation]] as <math>a\times10^{2n}</math> where <math> 1\leq a<100</math> and ''n'' is an integer, and the range of possible square roots is <math>\sqrt a\times10^n</math> where <math>1\leq \sqrt a<10</math>. 
====Scalar estimates====
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 (<math>\sqrt{10}\approx3.16</math>) times <math>10^n</math> 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 <math>\sqrt{S} = \sqrt{a}\times10^n</math> can be estimated as<ref group="Note">The factors two and six are used because they approximate the [[geometric mean]]s of the lowest and highest possible values with the given number of digits: <math>\sqrt{\sqrt{1} \cdot \sqrt{10}} = \sqrt[4]{10} \approx 1.78 \,</math> and <math>\sqrt{\sqrt{10} \cdot \sqrt{100}} = \sqrt[4]{1000} \approx 5.62 \,</math>.</ref>
<math display="block"> \sqrt{S} \approx \begin{cases}
2 \cdot 10^n & \text{if } a < 10, \\
6 \cdot 10^n & \text{if } a \geq 10.
\end{cases}</math>

This estimate has maximum absolute error of <math>4\cdot10^n</math> at a = 100, and maximum relative error of 100% at a = 1.

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

====Linear estimates====
A better estimate, and the standard method used, is a linear approximation to the function <math>y = x^2</math> over a small arc.  If, as above, powers of the base are factored out of the number <math>S</math> and the interval reduced to <math>[1,100]</math>, 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 <math>y=8.7x-10</math>. Reordering, <math>x=0.115y+1.15</math>. Rounding the coefficients for ease of computation,
<math display="block">\sqrt{S}\approx (a/10+1.2)\cdot 10^n</math>

That is the best estimate ''on average'' that can be achieved with a single piece linear approximation of the function  y=x<sup>2</sup> in the interval <math>[1,100]</math>.  It has a maximum absolute error of 1.2 at a=100, and maximum relative error of 30% at S=1 and 10.<ref group=Note>The unrounded estimate has maximum absolute error of 2.65 at 100 and maximum relative error of 26.5% at y=1, 10 and 100</ref> 

To divide by 10, subtract one from the exponent of <math>a</math>, 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 <math>[1,100]</math> is less than one significant digit of precision; the relative error is greater than 1/2<sup>2</sup>, 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''&nbsp;=&nbsp;1 to ''y''&nbsp;=&nbsp;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,{{radic|100|2}}] and [{{radic|100|2}},100].  For three intervals, the bounds are the cube roots of 100: [1,{{radic|100|3}}], [{{radic|100|3}},({{radic|100|3}})<sup>2</sup>], and [({{radic|100|3}})<sup>2</sup>,100], etc.  For two intervals, {{radic|100|2}} = 10, a very convenient number.  Tangent lines are easy to derive, and are located at x = {{sqrt|1*{{sqrt|10}}}} and ''x'' = {{sqrt|10*{{sqrt|10}}}}.  Their equations are: <math>y = 3.56x - 3.16</math>  and <math>y = 11.2x - 31.6</math>.  Inverting, the square roots are: <math>x = 0.28y + 0.89</math> and <math>x = .089y + 2.8</math>. Thus for <math>S = a \cdot 10^{2n}</math>:
<math display="block"> \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}</math>

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.

====Hyperbolic estimates====
In some cases, hyperbolic estimates may be efficacious, because a hyperbola is also a convex curve and may lie along an arc of ''y''&nbsp;=&nbsp;''x''<sup>2</sup> better than a line. Hyperbolic estimates are more computationally complex, because they necessarily require a floating division.  A near-optimal hyperbolic approximation to ''x''<sup>2</sup> on the interval <math>[1,100]</math> is y&nbsp;=&nbsp;190/(10−x)&thinsp;−&thinsp;20.  Transposing, the square root is ''x''&nbsp;=&nbsp;10&thinsp;−&thinsp;190/(y+20).  Thus for <math>S=a\cdot 10^{2n}</math>:
<math display="block">\sqrt{S}\approx \left(10-\frac{190}{a+20}\right)\cdot 10^n</math> 

The division need be accurate to only one decimal digit, because the estimate overall is only that accurate, and can be done mentally.  This hyperbolic estimate is better on average than scalar or linear estimates. It has maximum absolute error of 1.58 at ''a''&nbsp;=&nbsp;100 and maximum relative error at ''a''&nbsp;=&nbsp;10, where the estimate of 3.67 is 16.0% higher than the root of 3.16.  If instead one performed Newton-Raphson iterations beginning with an estimate of 10, it would take two iterations to get to 3.66, matching the hyperbolic estimate.  For a more typical case like 75, the hyperbolic estimate of 8.00 is only 7.6% low, and 5 Newton-Raphson iterations starting at 75 would be required to obtain a more accurate result.

====Arithmetic estimates====
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 &times; 5), and 36 is a perfect square (6 &times; 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 halfway between the squares.  So any number between 25 and halfway to 36, which is 30.5, estimate 5; any number greater than 30.5 up to 36, estimate 6.<ref group=Note>If the number is exactly half way between two squares, like 30.5, guess the higher number which is 6 in this case</ref>  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:

{| class="wikitable" style="text-align:center;"
|-
! {{nobold|{{mvar|a}}}}
!rowspan=2| nearest square
!rowspan=2| <math>k = \sqrt{a}</math> est.
|- style="height:1em;"
|rowspan=2| 1
|- style="height:1em;"
|rowspan=2| 1 (= 1<sup>2</sup>)
|rowspan=2 style="text-align: right;"| 1
|- style="height:1em;"
|rowspan=2| 2.5
|- style="height:1em;"
|rowspan=2| 4 (= 2<sup>2</sup>)
|rowspan=2 style="text-align: right;"| 2
|- style="height:1em;"
|rowspan=2| 6.5
|- style="height:1em;"
|rowspan=2| 9 (= 3<sup>2</sup>)
|rowspan=2 style="text-align: right;"| 3
|- style="height:1em;"
|rowspan=2| 12.5
|- style="height:1em;"
|rowspan=2| 16 (= 4<sup>2</sup>)
|rowspan=2 style="text-align: right;"| 4
|- style="height:1em;"
|rowspan=2| 20.5
|- style="height:1em;"
|rowspan=2| 25 (= 5<sup>2</sup>)
|rowspan=2 style="text-align: right;"| 5
|- style="height:1em;"
|rowspan=2| 30.5
|- style="height:1em;"
|rowspan=2| 36 (= 6<sup>2</sup>)
|rowspan=2 style="text-align: right;"| 6
|- style="height:1em;"
|rowspan=2| 42.5
|- style="height:1em;"
|rowspan=2| 49 (= 7<sup>2</sup>)
|rowspan=2 style="text-align: right;"| 7
|- style="height:1em;"
|rowspan=2| 56.5
|- style="height:1em;"
|rowspan=2| 64 (= 8<sup>2</sup>)
|rowspan=2 style="text-align: right;"| 8
|- style="height:1em;"
|rowspan=2| 72.5
|- style="height:1em;"
|rowspan=2| 81 (= 9<sup>2</sup>)
|rowspan=2 style="text-align: right;"| 9
|- style="height:1em;"
|rowspan=2| 90.5
|- style="height:1em;"
|rowspan=2| 100 (= 10<sup>2</sup>)
|rowspan=2 style="text-align: right;"| 10
|- style="height1em;"
|rowspan=2| 100
|- style="height:1em;"
|colspan=2|
|}

The final operation is to multiply the estimate {{mvar|k}} by the power of ten divided by 2, so for <math>S = a\cdot 10^{2n}</math>,
<math display="block">\sqrt{S}\approx k\cdot 10^n</math>

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 <math>k^2 \le a < (k+1)^2</math>, then <math>\sqrt{a}</math> is approximately k plus a fraction, the difference between {{mvar|a}} and ''k''<sup>2</sup> divided by the difference between the two squares:
<math display="block">\sqrt{a}\approx k + R</math> where <math display="block">R = \frac{(a-k^2)}{(k+1)^2-k^2}</math>

The final operation, as above, is to multiply the result by the power of ten divided by 2;
<math display="block">\sqrt{S} = \sqrt{a}\cdot 10^n \approx (k + R)\cdot 10^n</math>

{{mvar|k}} is a decimal digit and {{mvar|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. {{math|1=75 = 75 &times; 10<sup>2{{dot}}0</sup>}}, so {{mvar|a}} is 75 and {{mvar|n}} is 0. From the multiplication tables, the square root of the mantissa must be 8 point ''something'' because {{mvar|a}} is between 8&times;8 = 64 and 9&times;9 = 81, so {{mvar|k}} is 8; ''something'' is the decimal representation of {{mvar|R}}.  The fraction {{mvar|R}} is 75 &minus; ''k''<sup>2</sup> = 11, the numerator, and 81 &minus; ''k''<sup>2</sup> = 17, the denominator.  11/17 is a little less than 12/18 = 2/3 = .67, so guess .66 (it's okay to guess here, the error is very small). The final estimate is {{math|1=8 + .66 = 8.66}}. {{sqrt|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.

===Binary estimates===
When working in the [[binary numeral system]] (as computers do internally), by expressing <math>S</math> as <math>a\times2^{2n}</math> where <math> 0.1_2\leq a<10_2</math>, the square root <math>\sqrt{S} = \sqrt{a}\times2^n</math> can be estimated as 
<math display="block"> \sqrt{S} \approx (0.485 + 0.485 \cdot a) \cdot 2^n</math>

which is the least-squares regression line to 3 significant digit coefficients. <math>\sqrt{a}</math> has maximum absolute error of 0.0408 at <math>a</math>=2, and maximum relative error of 3.0% at <math>a</math>=1.  A computationally convenient rounded estimate (because the coefficients are powers of 2) is:
:<math> \sqrt{S} \approx (0.5 + 0.5 \cdot a) \cdot 2^n</math><ref group=Note>This is incidentally the equation of the tangent line to ''y''&nbsp;=&nbsp;''x''<sup>2</sup> at ''y''&nbsp;=&nbsp;1.</ref>
which has maximum absolute error of 0.086 at 2 and maximum relative error of 6.1% at {{math|1=''a'' = 0.5}} and {{math|1=''a'' = 2.0}}.

For <math> S = 125348 = 1\;1110\;1001\;1010\;0100_2 = 1.1110\;1001\;1010\;0100_2\times2^{16}\, </math>, the binary approximation gives <math> \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</math>. <math>\sqrt{125348}=354.0</math>, 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<sup>4</sup>, so the estimate is good to 4+ bits.

An estimate for <math>a</math> good to 8 bits can be obtained by table lookup on the high 8 bits of <math>a</math>, 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<sub>2</sub> representing 1.8515625<sub>10</sub>, the entry is 10101110<sub>2</sub> representing 1.359375<sub>10</sub>, the square root of 1.8515625<sub>10</sub> to 8 bit precision (2+ decimal digits).

==Heron's method==  
{{redirect|Heron's method|the formula used to find the area of a triangle|Heron's formula}}
The first explicit [[algorithm]] for approximating <math>\ \sqrt{ S ~}\ </math> is known as '''Heron's method''', after the first-century [[Ancient Greece|Greek]] mathematician [[Hero of Alexandria]] who described the method in his [[AD 60]] work ''[[Hero of Alexandria#Bibliography|Metrica]]''.{{sfn|Heath|1921}} This method is also called the '''Babylonian method''' (not to be confused with the [[#History|Babylonian method for approximating hypotenuses]]), although there is no evidence that the method was known to [[Babylonia|Babylonians]].

Given a positive real number <math>S</math>, let {{math|''x''<sub>0</sub> > 0}} be any positive [[#Initial estimate|initial estimate]].
Heron's method consists in iteratively computing 
<math display="block">x_{n+1} = \frac{1}{2} \left(x_n + \frac{S}{x_n}\right),</math>
until the desired accuracy is achieved.
The sequence <math>\ \bigl(\ x_0,\ x_1,\ x_2,\ x_3,\ \ldots\ \bigr)\ </math> defined by this equation converges to <math>\ \lim_{n \to \infty} x_n = \sqrt{ S ~} ~.</math> 

This is equivalent to using [[Newton's method]] to solve <math> x^2 - S = 0 </math>. This algorithm is [[quadratic convergence|quadratically convergent]]: the number of correct digits of <math>x_n</math> roughly doubles with each iteration.

===Derivation===
The basic idea is that if <math>\ x\ </math> is an overestimate to the square root of a non-negative real number <math>\ S\ </math> then <math>\ \tfrac{\ S\ }{ x }\ </math> will be an underestimate, and vice versa, so the average of these two numbers may reasonably be expected to provide a better approximation (though the formal proof of that assertion depends on the [[inequality of arithmetic and geometric means]] that shows this average is always an overestimate of the square root, as noted in the article on [[Square root#Geometric construction of the square root|square roots]], thus assuring convergence). 

More precisely, if <math>\ x\ </math> is our initial guess of <math>\ \sqrt{ S ~}\ </math> and <math>\ \varepsilon\ </math> is the error in our estimate such that <math>\ S = \left( x + \varepsilon \right)^2\ ,</math> then we can expand the binomial as: 
<math display="block">\ \bigl(\ x + \varepsilon\ \bigr)^2 = x^2 + 2 x \varepsilon + \varepsilon^2 </math> 
and solve for the error term
:<math>\varepsilon = \frac{\ S - x^2\ }{\ 2 x + \varepsilon\ } \approx \frac{\ S - x^2\ }{ 2 x }\ ,</math> if we suppose that <math>\ \varepsilon \ll x ~</math> 
Therefore, we can compensate for the error and update our old estimate as 
<math display="block">\ x + \varepsilon\ \approx\ x + \frac{\ S - x^2\ }{ 2 x }\ =\ \frac{\ S + x^2\ }{ 2 x }\ =\ \frac{\ \frac{ S }{\ x\ } + x\ }{ 2 }\ \equiv\ x_\mathsf{revised} ~.</math> 
Since the computed error was not exact, this is not the actual answer, but becomes our new guess to use in the next round of correction. The process of updating is iterated until desired accuracy is obtained.

This algorithm works equally well in the [[p-adic number|{{mvar|p}}-adic number]]s, but cannot be used to identify real square roots with {{mvar|p}}-adic square roots; one can, for example, construct a sequence of rational numbers by this method that converges to {{math|+3}} in the reals, but to {{math|−3}} in the 2-adics.

===Example===
To calculate <math>\sqrt{ S\,}</math> for <math>S = 125348</math> to seven significant figures, use the rough estimation method above to get
<math display="block">\begin{alignat}{5}
x_0 & = 6 \cdot 10^2 &&&&= 600\\[0.3em]
x_1 & = \frac 12 \left( x_0 + \frac S{x_0} \right) && = \frac 12 \left(600\phantom{.1} + \frac{125348}{600}\right) && = 404.457 \approx 400 \\[0.3em]
x_2 & = \frac 12 \left( x_1 + \frac S{x_1} \right) && = \frac 12 \left(400\phantom{.1} + \frac{125348}{400}\right) && = 356.685 \approx 360 \\[0.3em]
x_3 & = \frac 12 \left( x_2 + \frac S{x_2} \right) && = \frac 12 \left(360\phantom{.1} + \frac{125348}{360}\right) && = 354.094 \approx 354.1 \\[0.3em]
x_4 & = \frac 12 \left( x_3 + \frac S{x_3} \right) && = \frac 12 \left(354.1 + \frac{125348}{354.1}\right) && = 354.045199
\end{alignat}</math>

Therefore <math>\sqrt{\,125348\,} \approx 354.0452</math> to seven decimal places.  (The true value is 354.0451948551....)  Notice that early iterations only needed to be computed to 1, 2 or 4 places to produce an accurate final answer.

===Convergence===
[[File:Babylonian_method_graphs.svg|thumb|300px|link={{filepath:Babylonian_method_graphs.svg}}|Semilog graphs comparing the speed of convergence of Heron's method to find the square root of {{math|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 [{{filepath:Babylonian_method_graphs.svg}} the SVG file,] hover over a graph to display its points.]]
<!--
[[Image:Babylonian method.svg|300px|right|thumb|
Graph charting the use of the Babylonian method for approximating a square root of 100 (±10) using starting values
<span style="color:#008">''x''<sub>0</sub>&nbsp;=&nbsp;50</span>,
<span style="color:#F00">''x''<sub>0</sub>&nbsp;=&nbsp;1</span>,
and <span style="color:#080">''x''<sub>0</sub>&nbsp;=&nbsp;&minus;5</span>. Note that a positive starting value yields the positive root, and a negative starting value the negative root.]]
-->

Suppose that <math>\ x_0 > 0 ~~\mathsf{ and }~~ S > 0 ~.</math> Then for any natural number <math>\ n : x_n > 0 ~.</math> Let the [[relative error]] in <math>\ x_n\ </math> be defined by 
<math display="block">\ \varepsilon_n = \frac{~ x_n\ }{\ \sqrt{ S ~}\ } - 1 > -1\ </math>
and thus
<math display="block">\ x_n = \sqrt{ S ~} \cdot \left( 1 + \varepsilon_n \right) ~.</math>

Then it can be shown that
<math display="block">\ \varepsilon_{n+1} = \frac {\varepsilon_n^2}{2 (1 + \varepsilon_n)} \geq 0 ~.</math>

And thus that
<math display="block">\ \varepsilon_{n+2} \leq \min \left\{\ \frac{\ \varepsilon_{n+1}^2\ }{ 2 }, \frac{\ \varepsilon_{n+1}\ }{ 2 }\ \right\}\ </math>
and consequently that convergence is assured, and [[quadratic convergence|quadratic]].

====Worst case for convergence====
If using the rough estimate above with the Babylonian method, then the least accurate cases in ascending order are as follows:
<math display="block">\begin{align}
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{align}</math>

Thus in any case,
<math display="block">\begin{align}
 \varepsilon_1 & \leq 2^{-2}. \\
 \varepsilon_2 & < 2^{-5} < 10^{-1} ~. \\
 \varepsilon_3 & < 2^{-11} < 10^{-3} ~. \\
 \varepsilon_4 & < 2^{-23} < 10^{-6} ~. \\
 \varepsilon_5 & < 2^{-47} < 10^{-14} ~. \\
 \varepsilon_6 & < 2^{-95} < 10^{-28} ~. \\
 \varepsilon_7 & < 2^{-191} < 10^{-57} ~. \\
 \varepsilon_8 & < 2^{-383} < 10^{-115} ~.
\end{align} </math>

Rounding errors will slow the convergence. It is recommended to keep at least one extra digit beyond the desired accuracy of the <math>\ x_n\ </math> being calculated, to avoid significant round-off error.

==Bakhshali method==
This method for finding an approximation to a square root was described in an Ancient [[India|Indian]] manuscript, called the [[Bakhshali manuscript]]. It is algebraically equivalent to two iterations of Heron's method and thus quartically convergent, meaning that the number of correct digits of the approximation roughly quadruples with each iteration.{{sfn|Bailey|Borwein|2012}} The original presentation, using modern notation, is as follows: To calculate <math>\sqrt{S}</math>, let <math>x_0^2</math> be the initial approximation to <math>S</math>. Then, successively iterate as:
<math display="block">\begin{align}
a_n &= \frac{S - x_n^2}{2 x_n},\\
x_{n+1} &= x_n + a_n,\\
x_{n+2} &= x_{n+1} - \frac{a_n^2}{2 x_{n+1}}.
\end{align}</math>

The values <math>x_{n+1}</math> and <math>x_{n+2}</math> are exactly the same as those computed by Heron's method.  To see this, the second Heron's method step would compute
<math display=block>x_{n+2} = \frac{x_{n+1}^2 + S}{2x_{n+1}} = x_{n+1} + \frac{S - x_{n+1}^2}{2x_{n+1}}</math>
and we can use the definitions of <math>x_{n+1}</math> and <math>a_n</math> to rearrange the numerator into:
<math display=block>\begin{align}
S - x_{n+1}^2 &= S - (x_n + a_n)^2\\
          &= S - x_n^2 - 2x_na_n - a_n^2\\
          &= S - x_n^2 - (S - x_n^2) - a_n^2\\
          &= - a_n^2.
\end{align}</math>

This can be used to construct a rational approximation to the square root by beginning with an integer.  If <math>x_0 = N</math> is an integer chosen so <math>N^2</math> is close to <math>S</math>, and <math>d = S - N^2</math> is the difference whose absolute value is minimized, then the first iteration can be written as:
<math display="block">\sqrt{S} \approx N + \frac{d}{2N} - \frac{d^2}{8N^3 + 4Nd} = \frac{8N^4 + 8N^2 d + d^2}{8N^3 + 4Nd} = \frac{N^4 + 6N^2S + S^2}{4N^3 + 4NS} = \frac{N^2(N^2 + 6S) + S^2}{4N(N^2 + S)}.</math>

The Bakhshali method can be generalized to the computation of an arbitrary root, including fractional roots.{{sfn|Simply Curious|2018}}

One might think the second half of the Bakhshali method could be used as a simpler form of Heron's iteration and used repeatedly, e.g.
<!--Yes, I know it's more conventional in mathematics to write
  -\frac{x}{y}, but \frac{-x}{y} makes the minus sign more legible.-->
<math display="block">\begin{align}
a_{n+1} &= \frac{-a_n^2}{2 x_{n+1}}, &
x_{n+2} &= x_{n+1} + a_{n+1}, \\
a_{n+2} &= \frac{-a_{n+1}^2}{2 x_{n+2}}, &
x_{n+3} &= x_{n+2} + a_{n+2}, \text{ etc.}
\end{align}</math>
however, this is [[numerically unstable]].  Without any reference to the original input value <math>S</math>, the accuracy is limited by that of the original computation of <math>a_n</math>, and that rapidly becomes inadequate.

===Example===
Using the same example <math>S = 125348</math> as in the [[#Example|Heron's method example]], the first iteration gives
<!--This example is rounded after each step.  The first to 2 digits (multiples of 10), the second to 8 digits (5 decimals).  This causes some wiggle in the last digit compared to carrying higher precision between steps.-->
<math display="block">\begin{alignat}{3}
x_0 &= 600\\[1ex]
a_0 &= \frac{125348 - 600^2}{2 \times 600} &&= -195.5433 \approx -200\\[1ex]
x_1 &= 600 + (-200) &&= \phantom{-}400\\[1ex]
x_2 &= 400 - \frac{(-200)^2}{2 \times 400} &&= \phantom{-}350
\end{alignat}</math>

Likewise the second iteration gives
<math display="block">\begin{alignat}{3}
a_2 &= \frac{125348 - 350^2}{2 \times 350} &&=\phantom{00}4.06857\\[1ex]
x_3 &= 350 + 4.06857 &&= 354.06857\\[1ex]
x_4 &= 354.06857 - \frac{4.06857^2}{2 \times 354.06857} &&= 354.045194
\end{alignat}</math>
Unlike in Heron's method, <math>x_3</math> must be computed to 8 digits because the formula for <math>x_4</math> does not correct any error in <math>x_3</math>.

==Digit-by-digit calculation==
{{missing information|section|identity and origin of this algorithm|date=August 2023}}
This is a method to find each digit of the square root in a sequence. This method is based on the [[binomial theorem]] and basically an inverse algorithm solving <math>(x+y)^2 = x^2 + 2xy + y^2</math>. It is slower than the Babylonian method, but it has several advantages:
* It can be easier for manual calculations.
* Every digit of the root found is known to be correct, i.e., it does not have to be changed later.
* If the square root has an expansion that terminates, the algorithm terminates after the last digit is found. Thus, it can be used to check whether a given integer is a [[square number]].
* The algorithm works for any [[numeral system|base]], and naturally, the way it proceeds depends on the base chosen.
Disadvantages are:
* It becomes unmanageable for higher roots.
* It does not tolerate inaccurate guesses or sub-calculations; such errors lead to every following digit of the result being wrong, unlike with [[Newton's method]], which self-corrects any approximation errors. 
* While digit-by-digit calculation is efficient enough on paper, it is much too expensive for software implementations. Each iteration involves larger numbers, requiring more memory, but only advances the answer by one correct digit. Thus algorithm takes more time for each additional digit. 

[[Napier's bones]] include an aid for the execution of this algorithm. The shifting {{math|''n''}}th root algorithm is a generalization of this method.

===Basic principle===
First, consider the case of finding the square root of a number {{math|''S''}}, that is the square of a base-10 two-digit number {{math|''XY''}}, where {{math|''X''}} is the tens digit and {{math|''Y''}} is the units digit. Specifically:
<math display="block">S = \left(10X + Y\right)^2 = 100X^2 + 20XY + Y^2.</math>
{{math|''S''}} will consist of 3 or 4 decimal digits.

Now to start the digit-by-digit algorithm, we split the digits of {{math|''S''}} in two groups of two digits, starting from the right. This means that the first group will be of 1 or 2 digits. Then we determine the value of {{math|''X''}} as the largest digit such that {{math|''X''<sup>2</sup>}} is less than or equal to the first group.
We then compute the difference between the first group and 
{{math|''X''<sup>2</sup>}} and start the second iteration by concatenating the second group to it. This is equivalent to subtracting <math>100X^2</math> from {{math|''S''}}, and we're left with <math>S' = 20XY + Y^2</math>.
We divide {{math|''S'''}} by 10, then divide it by {{math|''2X''}} and keep the integer part to try and guess {{math|''Y''}}.
We concatenate {{math|''2X''}} with the tentative {{math|''Y''}} and multiply it by {{math|''Y''}}. If our guess is correct, this is equivalent to computing: <math display="block">(10(2X) + Y)Y = 20XY + Y^2 = S',</math> and so the remainder, that is the difference between {{math|''S'''}} and the result, is zero; if the result is higher than {{math|''S' ''}}, we lower our guess by 1 and try again until the remainder is 0.
Since this is a simple case where the answer is a perfect square root {{math|''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 {{math|''S''}} by expressing it as a sum of ''n'' positive numbers such that
<math display="block">S = \left(a_1 + a_2 + a_3 + \dots + a_n\right)^2.</math>

By repeatedly applying the basic identity 
<math display="block">(x+y)^2 = x^2 +2xy + y^2,</math> 
the right-hand-side term can be expanded as
<math display="block">\begin{align}
& (a_1+a_2+a_3+ \dotsb +a_n)^2 \\
=& \, a_1^2 + 2a_1a_2 + a_2^2 + 2(a_1+a_2) a_3 + a_3^2 + \dots + 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 + \dots + \left[2 \left(\sum_{i=1}^{n-1} a_i\right) + a_n\right] a_n.
\end{align}
</math>

This expression allows us to find the square root by sequentially guessing the values of <math>a_i</math>s. Suppose that the numbers <math>a_1, \ldots, a_{m-1}</math> have already been guessed, then the ''m''-th term of the right-hand-side of above summation is given by <math>Y_{m} = \left[2 P_{m-1} + a_{m}\right] a_{m},</math> where <math display="inline">P_{m-1} = \sum_{i=1}^{m-1} a_i</math> is the approximate square root found so far. Now each new guess <math>a_m</math> should satisfy the recursion 
<math display="block">X_{m} = X_{m-1} - Y_{m},</math> 
where <math>X_{m}</math> is the sum of all the terms after <math>Y_{m}</math>, i.e. the remainder, such that <math>X_m \geq 0</math> for all <math>1\leq m\leq n,</math> with initialization <math>X_0 = S.</math> When <math>X_n = 0,</math> the exact square root has been found; if not, then the sum of the <math>a_i</math>s gives a suitable approximation of the square root, with <math>X_n</math> being the approximation error.

For example, in the decimal number system we have 
<math display="block">S = \left(a_1 \cdot 10^{n-1} + a_2 \cdot 10^{n-2} + \cdots + a_{n-1} \cdot 10 + a_n \right)^2,</math> 
where <math>10^{n-i}</math> are place holders and the coefficients <math>a_i \in \{0,1,2,\ldots,9\}</math>. At any m-th stage of the square root calculation, the approximate root found so far, <math>P_{m-1}</math> and the summation term <math>Y_m</math> are given by 
<math display="block">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},</math>
<math display="block">Y_m = \left[2P_{m-1} + a_m \cdot 10^{n-m}\right] 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)}.</math>

Here since the place value of <math>Y_m</math> is an even power of 10, we only need to work with the pair of most significant digits of the remainder <math>X_{m-1}</math>, whose first term is <math>Y_{m}</math>, 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 <math>a_i</math> is searched from a smaller set of binary digits {0,1}. This makes the computation faster since at each stage the value of <math>Y_m</math> is either <math>Y_m = 0</math> for <math>a_m = 0</math> or <math>Y_m = 2 P_{m-1} + 1</math> for <math>a_m = 1</math>. The fact that we have only two possible options for <math>a_m</math> also makes the process of deciding the value of <math>a_m</math> at ''m''-th stage of calculation easier. This is because we only need to check if <math>Y_m \leq X_{m-1}</math> for <math>a_m = 1.</math> If this condition is satisfied, then we take <math>a_m = 1</math>; if not then <math>a_m = 0.</math> Also, the fact that multiplication by 2 is done by left bit-shifts helps in the computation.

===Decimal (base 10)===
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:

# 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'''''.
# 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 <math>x(20p + x) \le c</math>. We will use a new variable '''''y''''' = ''x''(20''p'' + ''x'').
#** Note: 20''p'' + ''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 <math>x</math> 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''.
# Subtract ''y'' from ''c'' to form a new remainder.
# 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.

====Examples====
'''Find the square root of 152.2756.'''

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

===Binary numeral system (base 2)===
This section uses the formalism from [[#Digit-by-digit calculation|the digit-by-digit calculation section above]], with the slight variation that we let <math>N^2 = (a_n + \dotsb + a_0)^2</math>, with each <math>a_m = 2^m</math> or <math>a_m = 0</math>.<br />
We iterate all <math>2^m</math>, from <math>2^n</math> down to <math>2^0</math>, and build up an approximate solution <math>P_m = a_n + a_{n-1} + \ldots + a_m</math>, the sum of all <math>a_i</math> for which we have determined the value.<br />
To determine if <math>a_m</math> equals <math>2^m</math> or <math>0</math>, we let <math>P_m = P_{m+1} + 2^m</math>. If <math>P_m^2 \leq N^2</math> (i.e. the square of our approximate solution including <math>2^m</math> does not exceed the target square) then <math>a_m = 2^m</math>, otherwise <math>a_m = 0</math> and <math>P_m = P_{m+1}</math>.<br />
To avoid squaring <math>P_m</math> in each step, we store the difference <math>X_m = N^2 - P_m^2</math> and incrementally update it by setting <math>X_m = X_{m+1} - Y_m</math> with <math>Y_m = P_m^2 - P_{m+1}^2 = 2P_{m+1}a_m + a_m^2</math>.<br />
Initially, we set <math>a_n = P_n = 2^n</math> for the largest <math>n</math> with <math>(2^n)^2 = 4^n \leq N^2</math>.

As an extra optimization, we store <math>P_{m+1}2^{m+1}</math> and <math>(2^m)^2</math>, the two terms of <math>Y_m</math> in case that <math>a_m</math> is nonzero, in separate variables <math>c_m</math>, <math>d_m</math>:
<math display="block">c_m = P_{m+1}2^{m+1}</math>
<math display="block">d_m = (2^m)^2</math>
<math display="block">
Y_m = \begin{cases}
c_m + d_m & \text{if } a_m = 2^m \\
0 & \text{if } a_m = 0
\end{cases}
</math>

<math>c_m</math> and <math>d_m</math> can be efficiently updated in each step:
<math display="block">c_{m-1} = P_m2^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}
</math>
<math display="block">d_{m-1} = \frac{d_m}{4}</math>

Note that:
<math display="block">c_{-1} = P_0 2^0 = P_0 = N,</math> which is the final result returned in the function below.

{{anchor|Woo Abacus}}An implementation of this algorithm in C:{{sfn|Guy|UKC|1985}}

<syntaxhighlight lang="c" line="1">
int32_t isqrt(int32_t n) {
    assert(("sqrt input should be non-negative", n > 0));

    // X_(n+1)
    int32_t x = n;

    // c_n
    int32_t c = 0;

    // d_n 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_(m+1) ≥ Y_m then a_m = 2^m
            x -= c + d;        // X_m = X_(m+1) - Y_m
            c = (c >> 1) + d;  // c_(m-1) = c_m/2 + d_m (a_m is 2^m)
        }
        else {
            c >>= 1;           // c_(m-1) = c_m/2      (aₘ is 0)
        }
        d >>= 2;               // d_(m-1) = d_m/4
    }
    return c;                  // c_(-1)
}
</syntaxhighlight>

Faster algorithms, in binary and decimal or any other base, can be realized by using lookup tables—in effect trading [[space–time tradeoff|more storage space for reduced run time]].{{sfn|Steinarson|Corbit|Hendry|2003}}

==Exponential identity==
[[calculator|Pocket calculator]]s typically implement good routines to compute the [[exponential function]] and the [[natural logarithm]], and then compute the square root of ''S'' using the identity found using the properties of logarithms (<math>\ln x^n = n \ln x</math>) and exponentials (<math>e^{\ln x} = x</math>):{{cn|date=May 2020}}
<math display="block">\sqrt{S} = e^{\frac{1}{2}\ln S}.</math>
The denominator in the fraction corresponds to the ''n''th root. In the case above the denominator is 2, hence the equation specifies that the square root is to be found. The same identity is used when computing square roots with [[logarithm table]]s or [[slide rule]]s.

==A two-variable iterative method==
This method is applicable for finding the square root of <math>0 < S < 3 \,\!</math> and converges best for <math>S \approx 1</math>.
This, however, is no real limitation for a computer-based calculation, as in base 2 floating-point and fixed-point representations, it is trivial to multiply <math>S \,\!</math> by an integer power of 4, and therefore <math>\sqrt{S}</math> by the corresponding power of 2, by changing the exponent or by shifting, respectively. Therefore, <math>S \,\!</math> can be moved to the range <math display="inline">\tfrac{1}{2} \le S <2</math>. Moreover, the following method does not employ general divisions, but only additions, subtractions, multiplications, and divisions by powers of two, which are again trivial to implement. A disadvantage of the method is that numerical errors accumulate, in contrast to single variable iterative methods such as the Babylonian one.

The initialization step of this method is
<math display="block">\begin{align} a_0 &= S \\ c_0 &= S - 1 \end{align}</math>
while the iterative steps read
<math display="block">\begin{align}
a_{n+1} &= a_n - a_n c_n / 2 \\
c_{n+1} &= c_n^2 (c_n - 3) / 4
\end{align}</math>
Then, <math>a_n \to \sqrt{S}</math> (while <math>c_n \to 0</math>).

The convergence of <math>c_n \,\!</math>, and therefore also of <math>a_n \,\!</math>, is quadratic.

The proof of the method is rather easy. First, rewrite the iterative definition of <math>c_n </math> as
<math display="block">1 + c_{n+1} = (1 + c_n) (1 - \tfrac 1 2 c_n)^2 .</math>
Then it is straightforward to prove by induction that
<math display="block">S (1 + c_n) = a_n^2</math>
and therefore the convergence of <math>a_n \,\!</math> to the desired result <math>\sqrt{S}</math> is ensured by the convergence of <math>c_n \,\!</math> to 0, which in turn follows from <math>-1 < c_0 < 2 \,\!</math>.

This method was developed around 1950 by [[Maurice Wilkes|M. V. Wilkes]], [[David Wheeler (computer scientist)|D. J. Wheeler]] and [[Stanley Gill|S. Gill]]{{sfn|Wilkes|Wheeler|Gill|1951}} for use on [[Electronic Delay Storage Automatic Calculator|EDSAC]], one of the first electronic computers.{{sfn|Campbell-Kelly|2009}} The method was later generalized, allowing the computation of non-square roots.{{sfn|Gower|1958}}

==Iterative methods for reciprocal square roots==
The following are iterative methods for finding the reciprocal square root of ''S'' which is <math>1/\sqrt{S}</math>. Once it has been found, find <math>\sqrt{S}</math> by simple multiplication: <math>\sqrt{S} = S \cdot (1/\sqrt{S})</math>. These iterations involve only multiplication, and not division. They are therefore faster than the [[Methods of computing square roots#Heron's method|Babylonian method]]. However, they are not stable. If the initial value is not close to the reciprocal square root, the iterations will diverge away from it rather than converge to it. It can therefore be advantageous to perform an iteration of the Babylonian method on a rough estimate before starting to apply these methods.

*Applying [[Newton's method]] to the equation <math>(1/x^2) - S = 0</math> produces a method that converges quadratically using three multiplications per step: <math display="block">x_{n+1} = \frac{x_n}{2} \cdot (3 - S \cdot x_n^2) = x_n \cdot \left(\frac{3}{2} - \frac{S}{2}\cdot x_n^2\right).</math>
*Another iteration is obtained by [[Halley's method]], which is the [[Householder's method]] of order two. This [[rate of convergence|converges cubically]], but involves five multiplications per iteration:{{cn|date=August 2019}} <math display="block">y_n = S \cdot x_n^2,</math> and <math display="block">x_{n+1} = \frac{x_n}{8} \cdot (15 - y_n \cdot (10 - 3 \cdot y_n)) =  x_n \cdot \left(\frac{15}{8} - y_n \cdot \left(\frac{10}{8} - \frac{3}{8} \cdot y_n\right)\right).</math>
* If doing [[fixed-point arithmetic]], the multiplication by 3 and division by 8 can implemented using shifts and adds.  If using floating-point, Halley's method can be reduced to four multiplications per iteration by precomputing <math display="inline">\sqrt{3/8} S</math> and adjusting all the other constants to compensate: <math display="block">y_n = \sqrt\frac{3}{8} S \cdot x_n^2,</math> and <math display="block">x_{n+1} = x_n \cdot \left(\frac{15}{8} - y_n \cdot \left(\sqrt\frac{25}{6} - y_n\right)\right).</math>
===Goldschmidt's algorithm===

Goldschmidt's algorithm is an extension of [[Goldschmidt division]], named after Robert Elliot Goldschmidt,<ref>{{cite thesis |first=Robert E. |last=Goldschmidt |title=Applications of Division by Convergence |series=M.Sc. dissertation |publisher=M.I.T. |date=1964 |oclc=34136725 |url=http://dspace.mit.edu/bitstream/handle/1721.1/11113/34136725-MIT.pdf |access-date=2015-09-15 |archive-date=2015-12-10 |archive-url=https://web.archive.org/web/20151210223340/http://dspace.mit.edu/bitstream/handle/1721.1/11113/34136725-MIT.pdf |url-status=live }}</ref><ref>{{cite journal |title=Authors |url=https://ieeexplore.ieee.org/document/5392026 |journal=IBM Journal of Research and Development | year=1967 | volume=11 | pages=125–127 | doi=10.1147/rd.111.0125 |archive-url=https://web.archive.org/web/20180718114413/https://ieeexplore.ieee.org/stamp/stamp.jsp?arnumber=5392026 |archive-date=18 July 2018}}</ref> which can be used to calculate square roots. Some computers use Goldschmidt's algorithm to simultaneously calculate <math>\sqrt{S}</math> and <math>1 / \sqrt{S}</math>.
Goldschmidt's algorithm finds <math>\sqrt{S}</math> 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.{{sfn|Markstein|2004}}

The first way of writing Goldschmidt's algorithm begins
: <math>b_0 = S</math>
: <math>Y_0 \approx 1/\sqrt{S}</math> (typically using a table lookup)
: <math>y_0 = Y_0</math>
: <math>x_0 = S y_0</math>
and iterates
<math display="block">\begin{align}
b_{n+1} &= b_n Y_n^2 \\
Y_{n+1} &= \tfrac 1 2 (3 - b_{n+1}) \\
x_{n+1} &= x_n Y_{n+1} \\
y_{n+1} &= y_n Y_{n+1}
\end{align}</math>
until <math>b_i</math> is sufficiently close to 1, or a fixed number of iterations.  The iterations converge to
<math display="block">\lim_{n \to \infty} x_n = \sqrt S,</math> and
<math display="block">\lim_{n \to \infty} y_n = 1/\sqrt S.</math>
Note that it is possible to omit either <math>x_n</math> and <math>y_n</math> from the computation, and if both are desired then <math>x_n = S y_n</math> may be used at the end rather than computing it through in each iteration.

A second form, using [[fused multiply-add]] operations, begins
: <math>y_0 \approx 1/\sqrt{S}</math> (typically using a table lookup)
: <math>x_0 = S y_0</math>
: <math>h_0 = \tfrac 1 2 y_0</math>
and iterates
<math display="block">\begin{align}
r_n &= 0.5 - x_n h_n \\
x_{n+1} &= x_n + x_n r_n \\
h_{n+1} &= h_n + h_n r_n
\end{align}</math>
until <math>r_i</math> is sufficiently close to 0, or a fixed number of iterations.  This converges to
<math display="block">\lim_{n \to \infty} x_n = \sqrt S,</math> and
<math display="block">\lim_{n \to \infty} 2h_n = 1/\sqrt S.</math>

==Taylor series==
If ''N'' is an approximation to <math>\sqrt{S}</math>, a better approximation can be found by using the [[Taylor series]] of the [[square root]] function:
<math display="block">\sqrt{N^2+d} = N\sum_{n=0}^\infty \frac{(-1)^{n}(2n)!}{(1-2n)n!^2 4^n}\frac{d^n}{N^{2n}} = N\left(1 + \frac{d}{2N^2} - \frac{d^2}{8N^4} + \frac{d^3}{16N^6} - \frac{5d^4}{128N^8} + \cdots\right)</math>

As an iterative method, the [[order of convergence]] is equal to the number of terms used. With two terms, it is identical to the [[#Heron's method|Babylonian method]]. With three terms, each iteration takes almost as many operations as the [[#Bakhshali method |Bakhshali approximation]], but converges more slowly.{{Citation needed|date=September 2017}} Therefore, this is not a particularly efficient way of calculation. To maximize the rate of convergence, choose {{mvar|N}} so that <math> \frac{|d|}{N^2} \,</math> is as small as possible.

==Continued fraction expansion==
{{see also|Solving quadratic equations with continued fractions}}
The [[continued fraction]] representation of a real number can be used instead of its decimal or binary expansion and this representation has the property that the square root of any rational number (which is not already a perfect square) has a periodic, repeating expansion, similar to how rational numbers have repeating expansions in the decimal notation system.

[[Quadratic irrational]]s (numbers of the form <math>\frac{a+\sqrt{b}}{c}</math>, where ''a'', ''b'' and ''c'' are integers), and in particular, square roots of integers, have [[periodic continued fraction]]s. Sometimes what is desired is finding not the numerical value of a square root, but rather its [[continued fraction]] expansion, and hence its rational approximation. Let ''S'' be the positive number for which we are required to find the square root. Then assuming ''a'' to be a number that serves as an initial guess and ''r'' to be the remainder term, we can write <math>S = a^2 + r. </math> Since we have <math>S - a^2 = (\sqrt{S} + a)(\sqrt{S} - a) = r</math>, we can express the square root of ''S'' as
<math display="block"> \sqrt{S} = a + \frac{r}{a + \sqrt{S}}. </math>

By applying this expression for <math>\sqrt{S}</math> to the denominator term of the fraction, we have
<math display="block"> \sqrt{S} = a + \frac{r}{a + (a + \frac{r}{a + \sqrt{S}})} = a + \frac{r}{2a + \frac{r}{a + \sqrt{S}}}. </math>
<div style="float:right; padding:1em; margin:0 0 0 1em; width:450px; border:1px solid; background:#faf3ee;">
{{center|'''<big>Compact notation</big>'''}}
The numerator/denominator expansion for continued fractions (see left) is cumbersome to write as well as to embed in  text formatting systems. So mathematicians have devised several alternative notations, like<ref>see: [[Generalized continued fraction#Notation]]</ref> 
<math display="block">
\sqrt{S} = a+
\frac{r}{2a+}\,
\frac{r}{2a+}\,
\frac{r}{2a+}\cdots
</math>

When <math>r = 1</math> throughout, an even more compact notation is:<ref>see: [[Continued fraction#Notations]]</ref>
<math display="block">[a; 2a, 2a, 2a, \cdots] </math>
For repeating continued fractions (which all square roots of non-perfect squares do), the repetend is represented only once, with an overline to signify a non-terminating repetition of the overlined part:<ref>see: [[Periodic continued fraction]]</ref> 
<math display="block">[a;\overline{2a}]</math>

For {{sqrt|2}}, the value of <math>a</math> is 1, so its representation is:
<math display="block">[1;\overline{2}]</math>
</div>

Proceeding this way, we get a [[generalized continued fraction]] for the square root as
<math display="block"> \sqrt{S} = a + \cfrac{r}{2a + \cfrac{r}{2a + \cfrac{r}{2a + \ddots}}}</math> 

The first step to evaluating such a fraction{{sfn|Sardina|2007|loc=2.3j on p.10}} to obtain a root is to do numerical substitutions for the root of the number desired, and number of denominators selected. For example, in canonical form, <math>r</math> is 1 and for {{sqrt|2}}, <math>a</math> is 1, so the numerical continued fraction for 3 denominators is: 
<math display="block"> \sqrt{2} \approx 1 + \cfrac{1}{2 + \cfrac{1}{2 + \cfrac{1}{2}}}</math>

Step 2 is to reduce the continued fraction from the bottom up, one denominator at a time, to yield a rational fraction whose numerator and denominator are integers. The reduction proceeds thus (taking the first three denominators):
<math display="block">\begin{align}
 1 + \cfrac{1}{2 + \cfrac{1}{2 + \cfrac{1}{2}}}
 &= 1 + \cfrac{1}{2 + \cfrac{1}{\frac{5}{2}}} \\
 &= 1 + \cfrac{1}{2 + \cfrac{2}{5}}
 = 1 + \cfrac{1}{\frac{12}{5}} \\
 &= 1 + \cfrac{5}{12}
 = \frac{17}{12}
\end{align}</math>

Finally (step 3), divide the numerator by the denominator of the rational fraction to obtain the approximate value of the root:
<math display="block">17 \div 12 = 1.42</math> rounded to three digits of precision.

The actual value of {{sqrt|2}} is 1.41 to three significant digits. The relative error is 0.17%, so the rational fraction is good to almost three digits of precision.  Taking more denominators gives successively better approximations: four denominators yields the fraction <math>\frac{41}{29} = 1.4137</math>, good to almost 4 digits of precision, etc.

The following are examples of square roots, their simple continued fractions, and their first terms {{mdash}} called ''convergents'' {{mdash}} up to and including denominator 99:
{| class="wikitable"
|-
! {{sqrt|S}}
! ~decimal
! continued fraction
! convergents
|-
| {{sqrt|2}}
| 1.41421
| style="text-align: left;" | <math>[1;\overline {2}]</math>
| <math>\frac{3}{2}, \frac{7}{5}, \frac{17}{12}, \frac{41}{29}, \frac{99}{70}</math>
|-
| {{sqrt|3}}
| 1.73205
| style="text-align: left;" | <math>[1;\overline {1,2}]</math>
| <math>\frac{2}{1}, \frac{5}{3}, \frac{7}{4}, \frac{19}{11}, \frac{26}{15}, \frac{71}{41}, \frac{97}{56}</math>
|-
| {{sqrt|5}}
| 2.23607
| style="text-align: left;" | <math>[2;\overline {4}]</math>
| <math>\frac{9}{4}, \frac{38}{17}, \frac{161}{72}</math>
|-
| {{sqrt|6}}
| 2.44949
| style="text-align: left;" | <math>[2;\overline {2,4}]</math>
| <math>\frac{5}{2}, \frac{22}{9}, \frac{49}{20}, \frac{218}{89}</math>
|-
| {{sqrt|10}}
| 3.16228
| style="text-align: left;" | <math>[3;\overline {6}]</math>
| <math>\frac{19}{6}, \frac{117}{37}</math>
|-
| <math>\sqrt{\pi}</math>
| 1.77245
| style="text-align: left;" | <math>[1;1,3,2,1,1,6...]</math>
| <math>\frac{2}{1}, \frac{7}{4}, \frac{16}{9}, \frac{23}{13}, \frac{39}{22}</math>
|-
| <math>\sqrt{e}</math>
| 1.64872
| style="text-align: left;" | <math>[1;1,1,1,5,1,1...]</math>
| <math>\frac{2}{1}, \frac{3}{2}, \frac{5}{3}, \frac{28}{17}, \frac{33}{20}, \frac{61}{37}</math>
|-
| <math>\sqrt{\phi}</math>
| 1.27202
| style="text-align: left;" | <math>[1;3,1,2,11,3,7...]</math>
| <math>\frac{4}{3}, \frac{5}{4}, \frac{14}{11}</math>
|}

In general, the larger the denominator of a rational fraction, the better the approximation. It can also be shown that truncating a continued fraction yields a rational fraction that is the best approximation to the root of any fraction with denominator less than or equal to the denominator of that fraction {{mdash}} e.g., no fraction with a denominator less than or equal to 70 is as good an approximation to {{sqrt|2}} as 99/70.

==Approximations that depend on the floating point representation==
<!--

    This section needs a diagram or something similar,
    to make it more accessible to non-techies.

-->
A number is represented in a [[floating point]] format as <math>m\times b^p</math> which is also called [[scientific notation]]. Its square root is <math>\sqrt{m}\times b^{p/2}</math> and similar formulae would apply for cube roots and logarithms. On the face of it, this is no improvement in simplicity, but suppose that only an approximation is required: then just <math>b^{p/2}</math> is good to an order of magnitude. Next, recognise that some powers, {{mvar|p}}, will be odd, thus for 3141.59 = 3.14159{{x10^|3}} rather than deal with fractional powers of the base, multiply the mantissa by the base and subtract one from the power to make it even. The adjusted representation will become the equivalent of 31.4159{{x10^|2}} so that the square root will be {{radic|31.4159}}{{x10^|1}}.

If the integer part of the adjusted mantissa is taken, there can only be the values 1 to 99, and that could be used as an index into a table of 99 pre-computed square roots to complete the estimate. A computer using base sixteen would require a larger table, but one using base two would require only three entries: the possible bits of the integer part of the adjusted mantissa are 01 (the power being even so there was no shift, remembering that a [[Normalized number|normalised]] floating point number always has a non-zero high-order digit) or if the power was odd, 10 or 11, these being the first ''two'' bits of the original mantissa. Thus, 6.25 = 110.01 in binary, normalised to 1.1001 &times; 2<sup>2</sup> an even power so the paired bits of the mantissa are 01, while .625 = 0.101 in binary normalises to 1.01 &times; 2<sup>−1</sup> an odd power so the adjustment is to 10.1 &times; 2<sup>−2</sup> and the paired bits are 10. Notice that the low order bit of the power is echoed in the high order bit of the pairwise mantissa. An even power has its low-order bit zero and the adjusted mantissa will start with 0, whereas for an odd power that bit is one and the adjusted mantissa will start with 1. Thus, when the power is halved, it is as if its low order bit is shifted out to become the first bit of the pairwise mantissa.

A table with only three entries could be enlarged by incorporating additional bits of the mantissa. However, with computers, rather than calculate an interpolation into a table, it is often better to find some simpler calculation giving equivalent results. Everything now depends on the exact details of the format of the representation, plus what operations are available to access and manipulate the parts of the number. For example, [[Fortran]] offers an <code>EXPONENT(x)</code> function to obtain the power. Effort expended in devising a good initial approximation is to be recouped by thereby avoiding the additional iterations of the refinement process that would have been needed for a poor approximation. Since these are few (one iteration requires a divide, an add, and a halving) the constraint is severe.

Many computers follow the [[IEEE floating-point standard|IEEE]] (or sufficiently similar) representation, and a very rapid approximation to the square root can be obtained for starting Newton's method. The technique that follows is based on the fact that the floating point format (in base two) approximates the base-2 logarithm. That is <math>\log_2(m\times 2^p) = p + \log_2(m)</math>

So for a 32-bit single precision floating point number in IEEE format (where notably, the power has a [[Exponent bias|bias]] of 127 added for the represented form) you can get the approximate logarithm by interpreting its binary representation as a 32-bit integer, scaling it by <math>2^{-23}</math>, and removing a bias of 127, i.e.
<math display="block">x_\text{int} \cdot 2^{-23} - 127 \approx \log_2(x).</math>

For example, 1.0 is represented by a [[hexadecimal]] number 0x3F800000, which would represent <math>1065353216 = 127 \cdot 2^{23}</math> if taken as an integer. Using the formula above you get <math>1065353216 \cdot 2^{-23} - 127 = 0</math>, as expected from <math>\log_2(1.0)</math>. In a similar fashion you get 0.5 from 1.5 (0x3FC00000).

[[Image:Log2approx.png]]

To get the square root, divide the logarithm by 2 and convert the value back. The following program demonstrates the idea. The exponent's lowest bit is intentionally allowed to propagate into the mantissa.  One way to justify the steps in this program is to assume <math>b</math> is the exponent bias and <math>n</math> is the number of explicitly stored bits in the mantissa and then show that
<math display="block">\left(\left( \tfrac 1 2 \left(x_\text{int} / 2^n - b\right) \right) + b\right) \cdot 2^n = \tfrac 1 2 \left(x_\text{int} - 2^n\right) + \left(\tfrac 1 2 \left(b + 1\right)\right) \cdot 2^n.</math>

<syntaxhighlight lang="c">
/* Assumes that float is in the IEEE 754 single precision floating point format */
#include <stdint.h>
float sqrt_approx(float z)
{
	union { float f; uint32_t i; } val = {z};	/* Convert type, preserving bit pattern */
	/*
	 * To justify the following code, prove that
	 *
	 * ((((val.i / 2^m) - b) / 2) + b) * 2^m = ((val.i - 2^m) / 2) + ((b + 1) / 2) * 2^m)
	 *
	 * where
	 *
	 * b = exponent bias
	 * m = number of mantissa bits
	 */
	val.i -= 1 << 23;	/* Subtract 2^m. */
	val.i >>= 1;		/* Divide by 2. */
	val.i += 1 << 29;	/* Add ((b + 1) / 2) * 2^m. */

	return val.f;		/* Interpret again as float */
}
</syntaxhighlight>

The three mathematical operations forming the core of the above function can be expressed in a single line.  An additional adjustment can be added to reduce the maximum relative error.  So, the three operations, not including the cast, can be rewritten as
<syntaxhighlight lang="c">
	val.i = (1 << 29) + (val.i >> 1) - (1 << 22) + a;
</syntaxhighlight>

where ''a'' is a bias for adjusting the approximation errors. For example, with ''a'' = 0 the results are accurate for even powers of 2 (e.g. 1.0), but for other numbers the results will be slightly too big (e.g. 1.5 for 2.0 instead of 1.414... with 6% error). With ''a'' = &minus;0x4B0D2, the maximum relative error is minimized to ±3.5%.

If the approximation is to be used for an initial guess for [[Newton's method]] to the equation <math>(1/x^2) - S = 0</math>, then the reciprocal form shown in the following section is preferred.

===Reciprocal of the square root===
{{Main|Fast inverse square root}}
A variant of the above routine is included below, which can be used to compute the [[Multiplicative inverse|reciprocal]] of the square root, i.e., <math>x^{-1/2}</math> instead, was written by Greg Walsh. The integer-shift approximation produced a relative error of less than 4%, and the error dropped further to 0.15% with one iteration of [[Newton's method]] on the following line.{{sfn|Lomont|2003}} In computer graphics it is a very efficient way to normalize a vector.

<syntaxhighlight lang="c">
float invSqrt(float x) {
    float xhalf = 0.5f * x;
    union {
        float x;
        int i;
    } u;
    u.x = x;
    u.i = 0x5f375a86 - (u.i >> 1);
    /* The next line can be repeated any number of times to increase accuracy */
    u.x = u.x * (1.5f - xhalf * u.x * u.x);
    return u.x;
}
</syntaxhighlight>

Some VLSI hardware implements inverse square root using a second degree polynomial estimation followed by a [[Division algorithm#Goldschmidt division|Goldschmidt iteration]].{{sfn|Piñeiro|Díaz Bruguera|2002}}

==Negative or complex square==
If ''S''&nbsp;<&nbsp;0, then its principal square root is
<math display="block">\sqrt {S} = \sqrt {\vert S \vert} \, \, i \,.</math>

If ''S''&nbsp;=&nbsp;''a''+''bi'' where ''a'' and ''b'' are real and ''b''&nbsp;≠&nbsp;0, then its principal square root is
<math display="block">\sqrt {S} = \sqrt{\frac{\vert S \vert + a}{2}} \, + \, \sgn (b) \sqrt{\frac{\vert S \vert - a}{2}} \, \, i \,.</math>

This can be verified by squaring the root.{{sfn|Abramowitz|Stegun|1964|loc=Section 3.7.26}}{{sfn|Cooke|2008}}
Here
<math display="block">\vert S \vert = \sqrt{a^2 + b^2}</math>

is the [[absolute value|modulus]] of ''S''. The principal square root of a [[complex number]] is defined to be the root with the non-negative real part.

== See also ==
* [[Alpha max plus beta min algorithm]]
* [[nth root algorithm]]

== Notes ==
{{Reflist|group=Note|30em}}

== References ==
{{Reflist}}

== Bibliography==
{{Refbegin|30em}}
*{{cite book
 | last1=Abramowitz
 | first1=Miltonn
 | last2=Stegun
 | first2=Irene A.
 | title=Handbook of mathematical functions with formulas, graphs, and mathematical tables
 | publisher=Courier Dover Publications
 | year=1964
 | isbn=978-0-486-61272-0
 | page=17
}}
*{{cite news
 | last1=Bailey
 | first1=David
 | last2=Borwein
 | first2=Jonathan
 | author-link2=Jonathan Borwein
 | title=Ancient Indian Square Roots: An Exercise in Forensic Paleo-Mathematics
 | journal=American Mathematical Monthly
 | volume = 119
 | issue=8
 | pages=646–657
 | date=2012
 | url=http://crd-legacy.lbl.gov/~dhbailey/dhbpapers/india-sqrt.pdf
 | access-date=2017-09-14
}}
*{{cite journal
 | last=Campbell-Kelly
 | first=Martin
 | author-link=Martin Campbell-Kelly
 | title=Origin of Computing
 | journal=Scientific American
 | volume=301
 | date=September 2009
 | issue=3
 | pages=62–69
 | doi=10.1038/scientificamerican0909-62
 | jstor=26001527
 | pmid=19708529
 | bibcode=2009SciAm.301c..62C
 }}
*{{cite book
 | last=Cooke
 | first=Roger 
 | title=Classical algebra: its nature, origins, and uses
 | publisher=John Wiley and Sons
 | year=2008
 | isbn=978-0-470-25952-8
 | page=59
}}
*{{cite journal
 | last1=Fowler
 | first1=David
 | last2=Robson
 | first2=Eleanor
 | title=Square Root Approximations in Old Babylonian Mathematics: YBC 7289 in Context
 | journal=Historia Mathematica
 | volume=25
 | date=1998
 | issue=4
 | page=376
 | doi=10.1006/hmat.1998.2209
 | url=https://www.sciencedirect.com/science/article/pii/S0315086098922091/pdf?md5=7c3e116fac2deef21b107b78f3561862&pid=1-s2.0-S0315086098922091-main.pdf
 | doi-access=free
}}
*{{cite journal
 | last=Gower
 | first=John C.
 | title=A Note on an Iterative Method for Root Extraction
 | journal=The Computer Journal
 | volume=1
 | issue=3
 | year=1958
 | pages=142–143
 | url=https://www.researchgate.net/publication/31300990
 | doi=10.1093/comjnl/1.3.142
 | doi-access=free
 }}
* {{cite journal
 | last=Jackson
 | first=Terence
 | date=2011-07-01
 | title=95.42 Irrational square roots of natural numbers — a geometrical approach
 | url=https://www.cambridge.org/core/journals/mathematical-gazette/article/abs/9542-irrational-square-roots-of-natural-numbers-a-geometrical-approach/6B9D8EBFDCC016013D303AA78973429F
 | journal=The Mathematical Gazette
 | language=en
 | volume=95
 | issue=533
 | pages=327–330
 | doi=10.1017/S0025557200003193
 | s2cid=123995083
 | issn=0025-5572}}
*{{cite web
 | last1=Guy
 | first1=Martin
 | last2=UKC
 |url=http://medialab.freaknet.org/martin/src/sqrt/sqrt.c |archive-url=https://web.archive.org/web/20120306040058/http://medialab.freaknet.org/martin/src/sqrt/sqrt.c |archive-date=2012-03-06 
 | title = Fast integer square root by Mr. Woo's abacus algorithm (archived)
 | date = 1985
}}
*{{cite book
 | last = Heath
 | first = Thomas
 | author-link = Thomas Little Heath
 | title = A History of Greek Mathematics, Vol. 2
 | publisher = Clarendon Press
 | year = 1921
 | location = Oxford
 | pages = [https://archive.org/details/ahistorygreekma00heatgoog/page/n340 323]–324
 | url = https://archive.org/details/ahistorygreekma00heatgoog
}}
*{{cite web
 | last=Lomont
 | first=Chris
 | title = Fast Inverse Square Root
 | url = http://www.lomont.org/Math/Papers/2003/InvSqrt.pdf
 | date = 2003
}}
*{{cite conference
 | last=Markstein
 | first=Peter
 | title=Software Division and Square Root Using Goldschmidt's Algorithms
 | date=November 2004
 | conference=6th Conference on Real Numbers and Computers
 | location=[[Dagstuhl]], Germany |url=http://www.informatik.uni-trier.de/Reports/TR-08-2004/rnc6_12_markstein.pdf
 | citeseerx=10.1.1.85.9648
 | conference-url=http://cca-net.de/rnc6/
}}
*{{cite journal
 | last1=Piñeiro
 | first1=José-Alejandro
 | last2=Díaz Bruguera
 | first2=Javier
 | title=High-Speed Double-Precision Computationof Reciprocal, Division, Square Root, and Inverse Square Root
 | journal=IEEE Transactions on Computers
 | volume=51
 | issue=12
 | date=December 2002
 | pages=1377–1388
 | url=https://www.researchgate.net/publication/31300990
 | doi=10.1109/TC.2002.1146704
}}
*{{cite web
 | last=Sardina
 | first=Manny
 | url = https://studylib.net/doc/7979641/general-method-for-extracting-roots-using--folded
 | title = General Method for Extracting Roots using (Folded) Continued Fractions
 | year= 2007
 | location = Surrey (UK)
}}
*{{cite web
 | author=Simply Curious
 | url=https://simplycurious.blog/2018/06/05/bucking-down-to-the-bakhshali-manuscript/
 | title=Bucking down to the Bakhshali manuscript
 | date=5 June 2018
 | access-date=2020-12-21
 | website=Simply Curious blog
}}
*{{cite web
 | author=Steinarson
 | first1=Arne
 | last2=Corbit
 | first2=Dann
 | last3=Hendry
 | first3=Mathew
 | url = http://atoms.alife.co.uk/sqrt/SquareRoot.java
 | title = Integer Square Root function
 | year= 2003
}}
*{{cite book
 | last1 = Wilkes
 | first1 = M.V.
 | author-link1=Maurice Wilkes
 | last2 = Wheeler
 | first2 = D.J.
 | author-link2=David Wheeler (computer scientist)
 | last3 = Gill
 | first3 = S.
 | author-link3=Stanley Gill
 | title = The Preparation of Programs for an Electronic Digital Computer
 | publisher = Addison-Wesley
 | location = Oxford
 | language = English
 | year = 1951
 | oclc = 475783493
 | pages = [https://archive.org/details/ahistorygreekma00heatgoog/page/n340 323]–324
 | url = https://archive.org/details/ahistorygreekma00heatgoog
}}
{{Refend}}

== External links ==
* {{MathWorld|title=Square root algorithms|urlname=SquareRootAlgorithms}}
* [https://web.archive.org/web/20220410042055/http://www.afjarvis.staff.shef.ac.uk/maths/jarvisspec02.pdf Square roots by subtraction]
* [http://www.andrijar.com/algorithms/algorithms.htm#qusr Integer Square Root Algorithm by Andrija Radović]
* [http://www.hparchive.com/Journals/HPJ-1977-05.pdf Personal Calculator Algorithms I : Square Roots (William E. Egbert), ''Hewlett-Packard Journal'' (May 1977) : page 22]
* [http://www.calculatorsquareroot.com Calculator to learn the square root]

[[Category:Computer arithmetic algorithms]]
[[Category:Root-finding algorithms]]