This article has multiple issues. Please helpimprove it or discuss these issues on thetalk page.(Learn how and when to remove these messages) This articlemay containoriginal research. Pleaseimprove it byverifying the claims made and addinginline citations. Statements consisting only of original research should be removed.(January 2012) (Learn how and when to remove this message) This articlemay be too technical for most readers to understand. Pleasehelp improve it tomake it understandable to non-experts, without removing the technical details.(September 2012) (Learn how and when to remove this message) (Learn how and when to remove this message)

Square root algorithms compute the non-negativesquare root${\displaystyle \sqrt{S}}$ of a positivereal number${\displaystyle S}$.Since all square roots ofnatural numbers, other than ofperfect squares, areirrational,^[1]square roots can usually only be computed to some finite precision: thesealgorithms typically construct a series of increasingly accurateapproximations.

Most square root computation methods are iterative: after choosing a suitableinitial estimate of${\displaystyle \sqrt{S}}$, aniterative refinement is performed until some termination criterion is met.One refinement scheme isHeron's method, a special case ofNewton's method.If division is much more costly than multiplication, it may be preferable to compute theinverse square root instead.

Other methods are available to compute the square rootdigit by digit, or usingTaylor series.Rational approximations of square roots may be calculated usingcontinued 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, aninteger 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).

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

Heron's method from first century Egypt was the first ascertainable algorithm for computing square root.^[3]

Modern analytic methods began to be developed after introduction of theArabic numeral system to western Europe in the early Renaissance.^[4]

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

Many iterative square root algorithms require an initialseed value. The seed must be a non-zero positive number; it should be between 1 and${\displaystyle S}$, 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${\displaystyle x_0=1}$ (or${\displaystyle S}$), then approximately${\displaystyle \frac{1}{2}|{\log }_{2}S|}$ 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${\displaystyle [x_0,S/x_0]}$). The estimate is a specific value of a functional approximation to${\displaystyle f(x)=\sqrt{x}}$ over the interval. Obtaining a better estimate involves either obtaining tighter bounds on the interval, or finding a better functional approximation to${\displaystyle f(x)}$. 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 andmantissa are usually treated separately, as the number would be expressed in scientific notation.

Typically the number${\displaystyle S}$ is expressed inscientific notation as${\displaystyle a\times {10}^{2n}}$ where andn is an integer, and the range of possible square roots is${\displaystyle \sqrt{a}\times {10}^n}$ where.

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 1]}

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

A better estimate, and the standard method used, is alinear approximation to the function${\displaystyle y=x^2}$ over a small arc. If, as above, powers of the base are factored out of the numberS and the interval reduced to [1, 100], asecant 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 estimateon average that can be achieved with a single piece linear approximation of the function${\displaystyle y=x^2}$ in the interval [1, 100]. It has a maximum absolute error of 1.2 ata=100, and maximum relative error of 30% atS=1 and 10.^{[Note 2]}

To divide by 10, subtract one from the exponent ofa, 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 [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 fromy = 1 toy = 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⁠${\displaystyle x=\sqrt{1\sqrt{10}}}$⁠ and⁠${\displaystyle x=\sqrt{10\sqrt{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}}$:

The maximum absolute errors occur at the high points of the intervals, ata=10 and 100, and are 0.54 and 1.7 respectively. The maximum relative errors are at the endpoints of the intervals, ata=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 ofy =x² better than a line. Hyperbolic estimates are more computationally complex, because they necessarily require a floating division. A near-optimal hyperbolic approximation tox² on the interval [1, 100] is⁠${\displaystyle y=190/(10-x)-20}$⁠. Transposing, the square root is⁠${\displaystyle x=10-190/(y+20)}$⁠. Thus for${\displaystyle S=a\cdot {10}^{2n}}$:$${\displaystyle \sqrt{S}\approx \left(10-\frac{190}{a+20}\right)\cdot {10}^n}$$

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 ata = 100 and maximum relative error ata = 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.

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 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.^{[Note 3]} 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
	1 (= 12)	1
2.5
	4 (= 22)	2
6.5
	9 (= 32)	3
12.5
	16 (= 42)	4
20.5
	25 (= 52)	5
30.5
	36 (= 62)	6
42.5
	49 (= 72)	7
56.5
	64 (= 82)	8
72.5
	81 (= 92)	9
90.5
	100 (= 102)	10
100

The final operation is to multiply the estimatek 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, then${\displaystyle \sqrt{a}}$ is approximately k plus a fraction, the difference betweena andk² 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 andR 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.

When working in thebinary numeral system (as computers do internally), by expressingS as${\displaystyle a\times 2^{2n}}$ where, the square root${\displaystyle \sqrt{S}=\sqrt{a}\times 2^n}$ can be estimated as$${\displaystyle \sqrt{S}\approx (0.485+0.485a)\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.5a)\cdot 2^n}$^{[Note 4]}

which has maximum absolute error of 0.086 at 2 and maximum relative error of 6.1% ata = 0.5 anda = 2.0.

For$${\displaystyle S=125348=1111010011010{0100}_2=1.111010011010{0100}_2\times 2^{16},}$$ the binary approximation gives$${\displaystyle \sqrt{S}\approx (0.5+0.5a)\cdot 2^8=1.011101001101{0010}_2\times 10000{0000}_2=1.456\times 256=\mathrm{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 fora good to 8 bits can be obtained by table lookup on the high 8 bits ofa, 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).

The first explicitalgorithm for approximating${\displaystyle \sqrt{S}}$ is known asHeron's method, after the first-centuryGreek mathematicianHero of Alexandria who described the method in hisAD 60 workMetrica.^[3] This method is also called theBabylonian method (not to be confused with theBabylonian method for approximating hypotenuses), although there is no evidence that the method was known toBabylonians.

Given a positive real number${\displaystyle S}$, letx₀ > 0 be any positiveinitial estimate.Heron's method consists in iteratively computing$${\displaystyle x_{n+1}=\frac{1}{2}\left(x_n+\frac{S}{x_n}\right),}$$until the desired accuracy is achieved.The sequence${\displaystyle (x_0,x_1,x_2,x_3,\dots )}$ defined by this equation converges to${\displaystyle \underset{n\to \mathrm{\infty }}{lim}{x}_n=\sqrt{S}.}$

This is equivalent to usingNewton's method to solve${\displaystyle x^2-S=0}$. This algorithm isquadratically convergent: the number of correct digits of${\displaystyle x_n}$ roughly doubles with each iteration.^[5]

The basic idea is that if${\displaystyle x}$ is an overestimate to the square root of a positive real number${\displaystyle S}$ then${\displaystyle \frac{S}{x}}$ will be an underestimate, and vice versa, so the average of these two numbers may reasonably be expected to provide a better approximation. (Theformal proof of that assertion depends on theinequality of arithmetic and geometric means that shows this average is always an overestimate of the square root, as noted in the article onsquare roots, thus assuring convergence.)

More precisely, if${\displaystyle x}$ is our initial guess of${\displaystyle \sqrt{S}}$ and${\displaystyle \epsilon }$ is the error in our estimate such that${\displaystyle S={\left(x+\epsilon \right)}^2,}$ then we can expand the binomial as:$${\displaystyle (x+\epsilon {)}^2=x^2+2x\epsilon +{\epsilon }^2}$$ and solve for the error term

${\displaystyle \epsilon =\frac{S-x^2}{2x+\epsilon }\approx \frac{S-x^2}{2x},}$ if we suppose that${\displaystyle \epsilon \ll x}$

Therefore, we can compensate for the error and update our old estimate as$${\displaystyle x+\epsilon \approx x+\frac{S-x^2}{2x}=\frac{S+x^2}{2x}=\frac{\frac{S}{x}+x}{2}\equiv x_{\mathsf{r}\mathsf{e}\mathsf{v}\mathsf{i}\mathsf{s}\mathsf{e}\mathsf{d}}.}$$ 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 thep-adic numbers, but cannot be used to identify real square roots withp-adic square roots; one can, for example, construct a sequence of rational numbers by this method that converges to+3 in the reals, but to−3 in the 2-adics.

fromdecimalimportDecimal,localcontext,getcontextNumberType=int|float|Decimaldefsqrt_Heron(s:NumberType,precision:int|None=None,guess:NumberType|None=None)->Decimal:"""    Compute sqrt(s) using the Heron-Newton method with arbitrary precision.    :param s: Non-negative number whose square root to compute.    :param precision: Number of significant digits.        Defaults to the current decimal context.        Minimum supported precision is 2.        (Precision = 1 is disallowed to avoid rounding anomalies.)    :param guess: Initial guess. Defaults to s / 2.    :return: Approximation of sqrt(s) rounded to the specified precision.    """ifs==0:returnDecimal(0)s=Decimal(s)ifs<0:raiseValueError("sqrt(s) is not defined for negative numbers.")ifprecisionisNone:precision=getcontext().prec# use current global context if not specified# Silently enforce minimum precisionifprecision<2:precision=2ifguessisNone:guess=Decimal(s/2)guard=25# temporary extra digits for internal stabilitymax_iter=10_000# Local context: isolate precision changeswithlocalcontext()asctx:ctx.prec=precision+guardguess=(guess+s/guess)/2for_inrange(max_iter):next_guess=(guess+s/guess)/2# Stop when improvement is small enoughifguess-next_guess<Decimal(f"1e-{precision}"):breakguess=next_guesselse:raiseArithmeticError(f"Heron method did not converge within{max_iter} iterations")# Round to target precision (getting rid of guard)ctx.prec=precisionreturn+next_guess

The following example demonstrates the execution of thesqrt_Heron function with various inputs

print(f"1){sqrt_Heron(125348,precision=7,guess=600)}")print(f"2){sqrt_Heron(Decimal('3.1415926535897932384626433832795028841971693993'))}")print(f"3){sqrt_Heron(2,1_000_157)}")print(f"4){sqrt_Heron(2,10_000_005,1.414)}")print(f"5){sqrt_Heron(2,100_000_000,1)}")

This produces the following output:

1) 354.04522) 1.7724538509055160272981674833) 1.4142135623730950488016887242 ... 2697320257318491414938800048567428924) 1.4142135623730950488016887242 ... ... 8724805080541235727278721315897142625) 1.4142135623730950488016887242 ... ... ... 023678977744844723443287604232894971

Ad 1)

The computation of${\displaystyle \sqrt{s}}$ for${\displaystyle s=125348}$ to seven significant figures takes the following course:

Therefore${\displaystyle \sqrt{125348}\approx 354.0452}$ to seven significant figures (rounded).

Ad 2) Computation (in 6 iteration steps) of${\displaystyle \sqrt{\lfloor \pi \times {10}^{46}\rfloor \times {10}^{-46}}}$ to default precision.^{[Note 5]}

Ad 3) Computation (in 22 iteration steps) of${\displaystyle \sqrt{2}}$ to 1,000,157 digits.^[6]

Ad 4) Computation (in 23 iteration steps) of${\displaystyle \sqrt{2}}$ to 10,000,005 digits.^[7]

Ad 5) Computation (in 28 iteration steps) of${\displaystyle \sqrt{2}}$ to 100 million digits.^{[citation needed]}

It appears that for reasonable initial guesses not many iterations are needed.

Heron's method has the following property:

In plain words:Once the iteration produces a value greater than${\displaystyle \sqrt{S}}$ (which happens immediately if${\displaystyle x_0>\sqrt{S}}$, or after one step if), every following estimate remains above${\displaystyle \sqrt{S}}$ but gets smaller each time — so the sequence “slides down” toward${\displaystyle \sqrt{S}}$ and converges.

In line 41 of the program,guess is set to a value${\displaystyle \ge \sqrt{S}}$. Then in line 46 of the code,${\displaystyle {\delta }_n=x_n-x_{n+1}}$ cannot be negative.

By using the difference between successive estimates,

${\displaystyle {\delta }_n=x_n-x_{n+1}}$,

as the stopping criterion, the method ensures that the sequence of approximations${\displaystyle x_n}$ is converging toward the true value${\displaystyle \sqrt{S}}$. When the successive differences${\displaystyle {\delta }_n}$ become sufficiently small, the given goal is reached. The key insight is that theabsolute error

${\displaystyle {\epsilon }_n=x_n-\sqrt{S}}$,

is directly related to the size of the successive improvement${\displaystyle {\delta }_n}$. Specifically, for iterative methods that converge linearly or quadratically, there exists a constant such that

${\displaystyle {\epsilon }_{n+1}=C\cdot {\delta }_n}$.

This relationship implies that as${\displaystyle {\delta }_n}$ decreases, the absolute error${\displaystyle {\epsilon }_n}$ also becomes smaller. Therefore, stopping the iteration when${\displaystyle {\delta }_n}$ falls below a given threshold ensures that the actual error is within at most that threshold.

Suppose that${\displaystyle x_0>0\mathsf{a}\mathsf{n}\mathsf{d}S>0.}$ Then for any natural number${\displaystyle n:x_n>0.}$ Let therelative error in${\displaystyle x_n}$ be defined by$${\displaystyle {\epsilon }_n=\frac{x_n}{\sqrt{S}}-1>-1}$$and thus$${\displaystyle x_n=\sqrt{S}\cdot \left(1+{\epsilon }_n\right).}$$

Then it can be shown that$${\displaystyle {\epsilon }_{n+1}=\frac{{\epsilon }_n^2}{2(1+{\epsilon }_n)}\ge 0.}$$

And thus that$${\displaystyle {\epsilon }_{n+2}\le min\left\{\frac{{\epsilon }_{n+1}^2}{2},\frac{{\epsilon }_{n+1}}{2}\right\}}$$and consequently that convergence is assured, andquadratic.

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

Thus in any case,

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

When in the program above the estimates — the highlighted lines 41 and 43 —

guess=(guess+s/guess)/2# ...next_guess=(guess+s/guess)/2

are replaced by

guess*=(guess*guess+3*s)/(3*guess*guess+s)# ...next_guess=guess*(guess*guess+3*s)/(3*guess*guess+s)

the functionsqrt_Heron is transformed into an implementation ofHalley's method, where

${\displaystyle x_{n+1}=x_n\cdot \frac{x_n^2+3S}{3x_n^2+S}}$

As in Halley's method the estimates do not always move in one direction, line 46 will become

ifabs(guess-next_guess)<Decimal(f"1e-{precision}"):

Halley's method converges faster —therate of convergence to the root is cubic, which is better thanquadratic— iteration for iteration, but involves five multiplications per iteration (counting division as three multiplications). The five example computations are completed in 4, respectively 4, 14, 15, and 19 iteration steps. By contrast Heron's method only needs one division i.e. three multiplications, so Heron's method is slightly better in the long run.

This method for finding an approximation to a square root was described in an AncientIndian manuscript, called theBakhshali 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.^[8] The original presentation, using modern notation, is as follows: To calculate${\displaystyle \sqrt{S}}$, let${\displaystyle x_0^2}$ be the initial approximation to${\displaystyle S}$. Then, successively iterate as:

The values${\displaystyle x_{n+1}}$ and${\displaystyle x_{n+2}}$ are exactly the same as those computed by Heron's method. To see this, the second Heron's method step would compute$${\displaystyle x_{n+2}=\frac{x_{n+1}^2+S}{2x_{n+1}}=x_{n+1}+\frac{S-x_{n+1}^2}{2x_{n+1}}}$$and we can use the definitions of${\displaystyle x_{n+1}}$ and${\displaystyle a_n}$ to rearrange the numerator into:

This can be used to construct a rational approximation to the square root by beginning with an integer. If${\displaystyle x_0=N}$ is an integer chosen so${\displaystyle N^2}$ is close to${\displaystyle S}$, and${\displaystyle d=S-N^2}$ is the difference whose absolute value is minimized, then the first iteration can be written as:$${\displaystyle \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^{2}S+S^2}{4N^3+4NS}=\frac{N^2(N^2+6S)+S^2}{4N(N^2+S)}.}$$

The Bakhshali method can be generalized to the computation of an arbitrary root, including fractional roots.^[9]

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.however, this isnumerically unstable. Without any reference to the original input value${\displaystyle S}$, the accuracy is limited by that of the original computation of${\displaystyle a_n}$, and that rapidly becomes inadequate.

Using the same example${\displaystyle S=125348}$ as in theHeron's method example, the first iteration gives

Likewise the second iteration givesUnlike in Heron's method,${\displaystyle x_3}$ must be computed to 8 digits because the formula for${\displaystyle x_4}$ does not correct any error in${\displaystyle x_3}$.

This is a method to find each digit of the square root in a sequence. The technique comes from the work ofFrançois Viète, published c. 1600.^[10] This method is based on thebinomial theorem and is essentially an inverse algorithm solving${\displaystyle (x+y{)}^2=x^2+2xy+y^2}$. 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 asquare number.
• The algorithm works for anybase, 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 withNewton'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 shiftingnth root algorithm is a generalization of this method.

First, consider the case of finding the square root of a numberS, that is the square of a base-10 two-digit numberXY, whereX is the tens digit andY is the units digit. Specifically:$${\displaystyle S={\left(10X+Y\right)}^2=100X^2+20XY+Y^2.}$$S will consist of 3 or 4 decimal digits.

Now to start the digit-by-digit algorithm, we split the digits ofS 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 ofX as the largest digit such thatX² is less than or equal to the first group.We then compute the difference between the first group andX² and start the second iteration by concatenating the second group to it. This is equivalent to subtracting${\displaystyle 100X^2}$ fromS, and we're left with${\displaystyle S^{\prime }=20XY+Y^2}$.We divideS' by 10, then divide it by2X and keep the integer part to try and guessY.We concatenate2X with the tentativeY and multiply it byY. If our guess is correct, this is equivalent to computing:$${\displaystyle (10(2X)+Y)Y=20XY+Y^2=S^{\prime },}$$ and so the remainder, that is the difference betweenS' and the result, is zero; if the result is higher thanS', 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 rootXY, 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 ofS by expressing it as a sum ofn positive numbers such that$${\displaystyle S={\left(a_1+a_2+a_3+\cdots +a_n\right)}^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

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,\dots ,a_{m-1}}$ have already been guessed, then them-th term of the right-hand-side of the above summation is given by${\displaystyle Y_m=\left[2P_{m-1}+a_m\right]a_m,}$ where$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,}$$ where${\displaystyle X_m}$ is the sum of all the terms after${\displaystyle Y_m}$, i.e. the remainder, such that${\displaystyle X_m\ge 0}$ for all${\displaystyle 1\le m\le n,}$ with initialization${\displaystyle X_0=S.}$ When${\displaystyle X_n=0,}$ the exact square root has been found; if not, then the sum of the${\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 S={\left(a_1\cdot {10}^{n-1}+a_2\cdot {10}^{n-2}+\cdots +a_{n-1}\cdot 10+a_n\right)}^2,}$$ where${\displaystyle {10}^{n-i}}$ are place holders and the coefficients${\displaystyle a_i\in \{0,1,2,\dots ,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=\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)}.}$$

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 remainder${\displaystyle X_{m-1}}$, whose first term is${\displaystyle Y_m}$, 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}$ atm-th stage of calculation easier. This is because we only need to check if${\displaystyle Y_m\le 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 thelong 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 thecurrent valuec.
2. Findp,y andx, as follows:
   • Letp be thepart of the root found so far, ignoring any decimal point. (For the first step,p = 0.)
   • Determine the greatest digitx such that${\displaystyle x(20p+x)\le c}$. We will use a new variabley =x(20p +x).
     • Note: 20p +x is simply twicep, with the digitx appended to the right.
     • Note:x can be found by guessing whatc/(20·p) is and doing a trial calculation ofy, then adjustingx 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 nextp will be the oldp times 10 plusx.
3. Subtracty fromc 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 00Algorithm terminates: Answer=12.34

This section uses the formalism fromthe digit-by-digit calculation section above, with the slight variation that we let${\displaystyle N^2=(a_n+\cdots +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}+\dots +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\le 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\le 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 c_m}$ and${\displaystyle d_m}$ can be efficiently updated in each step:$${\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.

ThePython program computes${\displaystyle \mathrm{isqrt}(n)=\lfloor \sqrt{n}\rfloor .}$. The algorithm is adigit-by-digit (bit-by-bit) method forinteger square roots.^[11]

defisqrt(x:int)->int:assertx>=0,"sqrt input should be non-negative"op:int=x# X_(n+1)res:int=0# c_n# d_n which starts at the highest power of four <= none:int=1whileone<=op:one<<=2# Now 'one' is the largest power of four <= xone>>=2# for dₙ … d₀whileone!=0:ifop>=res+one:# if X_(m+1) ≥ Y_m then a_m = 2^mop-=res+one# X_m = X_(m+1) - Y_mres+=2*one# c_m = c_m + 2*d_mres//=2# c_(m-1) = c_m / 2one//=4# d_(m-1) = d_m / 4# c_(-1)returnres

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

This sectiondoes notcite anysources. Please helpimprove this section byadding citations to reliable sources. Unsourced material may be challenged andremoved.(May 2020) (Learn how and when to remove this message)

Pocket calculators typically implement good routines to compute theexponential function and thenatural logarithm, and then compute the square root ofS using the identity found using the properties of logarithms (${\displaystyle \ln x^n=n\ln x}$) and exponentials(${\displaystyle e^{\ln x}=x}$):$${\displaystyle \sqrt{S}=e^{\frac{1}{2}\ln S}.}$$The denominator in the fraction corresponds to thenth 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 withlogarithm tables orslide rules.

This method is applicable for finding the square root of and converges best for${\displaystyle S\approx 1}$.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${\displaystyle S}$ by an integer power of 4, and therefore${\displaystyle \sqrt{S}}$ by the corresponding power of 2, by changing the exponent or by shifting, respectively. Therefore,${\displaystyle S}$ can be moved to the range. 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 iswhile the iterative steps readThen,${\displaystyle a_n\to \sqrt{S}}$ (while${\displaystyle c_n\to 0}$).

The convergence of${\displaystyle c_n}$, and therefore also of${\displaystyle a_n}$, is quadratic.

The proof of the method is rather easy. First, rewrite the iterative definition of${\displaystyle c_n}$ as$${\displaystyle 1+c_{n+1}=(1+c_n)(1-\frac{1}{2}{c}_n{)}^2.}$$Then it is straightforward to prove by induction that$${\displaystyle S(1+c_n)=a_n^2}$$and therefore the convergence of${\displaystyle a_n}$ to the desired result${\displaystyle \sqrt{S}}$ is ensured by the convergence of${\displaystyle c_n}$ to 0, which in turn follows from.

This method was developed around 1950 byM. V. Wilkes,D. J. Wheeler andS. Gill^[13] for use onEDSAC, one of the first electronic computers.^[14] The method was later generalized, allowing the computation of non-square roots.^[15]

The following are iterative methods for finding the reciprocal square root ofS which is${\displaystyle 1/\sqrt{S}}$. Once it has been found, find${\displaystyle \sqrt{S}}$ by simple multiplication:${\displaystyle \sqrt{S}=S\cdot (1/\sqrt{S})}$. These iterations involve only multiplication, and not division. They are therefore faster than theBabylonian 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.