Square root algorithms compute the non-negative square root of a positive real number. Since all square roots of natural numbers, other than of perfect squaress, are irrational, square roots can usually only be computed to some finite precision: these algorithms typically construct a series of increasingly accurate approximations.

Most square root computation methods are iterative: after choosing a suitable initial estimate, an iterative refinement is performed until some termination criterion is met. One refinement scheme is 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 inverse square root instead.

Other methods are available to compute the square root digit by digit, or using Taylor series. Rational approximations of square roots may be calculated using 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 mathematicians calculated the square root of 2 to three sexagesimal

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

Modern analytic methods began to be developed after introduction of the Arabic numeral system to western Europe in the early Renaissance.

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

. 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 are usually treated separately, as the number would be expressed in scientific notation.

Decimal estimates

is expressed in scientific notation and n.

Scalar estimates

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.

can be estimated as

.

Linear estimates

A better estimate, and the standard method used, is a linear approximation over a small arc. If, as above, powers of the base are factored out of the number and the interval reduced to , 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.

That is the best estimate on average in the interval . It has a maximum absolute error of 1.2 at a=100, and maximum relative error of 30% at S=1 and 10.

To divide by 10, subtract one from the exponent of , 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 is less than one significant digit of precision; the relative error is greater than 1/22, 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 to 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. and . For three intervals, the bounds are the cube roots of 100: , and , etc. For two intervals, , a very convenient number. Tangent lines are easy to derive, and are located at and

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 better than a line. Hyperbolic estimates are more computationally complex, because they necessarily require a floating division. A near-optimal hyperbolic approximation to x2 on the interval is . Transposing, the square root is

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 and maximum relative error at , 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 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. 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:The final operation is to multiply the estimate

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

is approximately k plus a fraction, the difference between and k2 divided by the difference between the two squares:

is a decimal digit and 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.

Binary estimates

When working in the binary numeral system (as computers do internally), by expressing

. A computationally convenient rounded estimate (because the coefficients are powers of 2) is:

which has maximum absolute error of 0.086 at 2 and maximum relative error of 6.1% at and .

, so the estimate has an absolute error of 19 and relative error of 5.3%. The relative error is a little less than 1/24, so the estimate is good to 4+ bits.

An estimate for good to 8 bits can be obtained by table lookup on the high 8 bits of , 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 111011012 representing 1.851562510, the entry is 101011102 representing 1.35937510, the square root of 1.851562510 to 8 bit precision (2+ decimal digits).

Heron's method

The first explicit algorithm is known as Heron's method, after the first-century Greek mathematician Hero of Alexandria who described the method in his AD 60 work Metrica. This method is also called the Babylonian method (not to be confused with the Babylonian method for approximating hypotenuses), although there is no evidence that the method was known to Babylonians.

, let be any positive initial estimate

This is equivalent to using Newton's method. This algorithm is quadratically convergent roughly doubles with each iteration.

Derivation

will be an underestimate, and vice versa, so the average of these two numbers may reasonably be expected to provide a better approximation. (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 roots, thus assuring convergence.)

and solve for the error term

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 numberss, but cannot be used to identify real square roots with -adic square roots; one can, for example, construct a sequence of rational numbers by this method that converges to in the reals, but to in the 2-adics.

Implementation in Python (programming language)Python

Example computation

The following example demonstrates the execution of the sqrt_Heron function with various inputs This produces the following output: 1) 354.0452 2) 1.772453850905516027298167483 3) 1.4142135623730950488016887242 ... 269732025731849141493880004856742892 4) 1.4142135623730950488016887242 ... ... 872480508054123572727872131589714262 5) 1.4142135623730950488016887242 ... ... ... 023678977744844723443287604232894971

Ad 1)

to seven significant figures

to seven significant figures (rounded).

Ad 2) to default precision.

Ad 3) to 1,000,157 digits.

Ad 4) to 10,000,005 digits.

Ad 5) to 100 million digits.

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

Notes

Explanation of lines 41 and 46

Heron's method has the following property:

and converges.

In line 41 of the program, guess cannot be negative.

Justification of the stopping criterion

By using the difference between successive estimates,

1. ,

become sufficiently small, the given goal is reached. The key insight is that the absolute error

1. ,

such that

1. .

falls below a given threshold ensures that the actual error is within at most that threshold.

Convergence

Let the relative error

and consequently that convergence is assured, and quadratic.

Worst case for convergence

being calculated, to avoid significant round-off error.

Halley's method

When in the program above the estimates the highlighted lines 41 and 43 are replaced by the function sqrt_Heron is transformed into an implementation of Halley's method, where

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

Halley's method converges faster the rate of convergence to the root is cubic, which is better than quadratic 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.

Bakhshali method

This method for finding an approximation to a square root was described in an Ancient 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.

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

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 is numerically unstable, and that rapidly becomes inadequate.

Example

as in the Heron's method example, the first iteration gives

.

Digit-by-digit calculation

This technique comes from the work of François Viète, published c. 1600., and is based on the binomial theorem. 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 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 th root algorithm is a generalization of this method.

Basic principle

First, consider the case of finding the square root of a number , that is the square of a base-10 two-digit number , where is the tens digit and will consist of 3 or 4 decimal digits.

Now to start the digit-by-digit algorithm, we split the digits of 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 as the largest digit such that is less than or equal to the first group. We then compute the difference between the first group and from . We divide by 10, then divide it by and keep the integer part to try and guess . We concatenate with the tentative and multiply it by and so the remainder, that is the difference between and the result, is zero; if the result is higher than , 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 , 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 by expressing it as a sum of n

have already been guessed, then the m being the approximation error.

, at any m-th stage. The section below codifies this procedure.

at m 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:

1. Starting on the left, bring down the most significant (leftmost) pair of digits not yet used (if all the digits have been used, write "00") and write them to the right of the remainder from the previous step (on the first step, there will be no remainder). In other words, multiply the remainder by 100 and add the two digits. This will be the current value c.

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. We will use a new variable y = x(20p + x).

• Note: 20p + x is simply twice p, with the digit x appended to the right.

• Note: x can be found by guessing what c/(20·p) is and doing a trial calculation of y, then adjusting x upward or downward as necessary.

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

          1  2. 3  4 
       /
     \/  01 52.27 56

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

Binary numeral system (base 2)

This section uses the formalism from the digit-by-digit calculation section above. for which we have determined the value....

which is the final result returned in the function below.

Implementation

The Python. The algorithm is a digit-by-digit (bit-by-bit) method for integer square roots.

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

Exponential identity

Pocket calculatorss typically implement good routines to compute the exponential function and the natural logarithm, and then compute the square root of S) and exponentials The denominator in the fraction corresponds to the nth 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 tables or slide rules.

A two-variable iterative method

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

).

, is quadratic.

.

This method was developed around 1950 by M. V. Wilkes, D. J. Wheeler and S. Gill for use on EDSAC, one of the first electronic computers. The method was later generalized, allowing the computation of non-square roots.

Iterative methods for reciprocal square roots

The following are iterative methods for finding the reciprocal square root of S. These iterations involve only multiplication, and not division. They are therefore faster than the Babylonian method. However, they are not stable. If the initial value is not close to the 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

• Another iteration is obtained by Halley's method, which is the Householder's method of order two. This converges cubically, but involves five multiplications per iteration:

• If doing fixed-point arithmetic

Goldschmidt's algorithm

Goldschmidt's algorithm is an extension of Goldschmidt division, named after Robert Elliot Goldschmidt, 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.

The first way of writing Goldschmidt's algorithm begins

1. (typically using a table lookup)

may be used at the end rather than computing it through in each iteration.

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

1. (typically using a table lookup)

Taylor series

If N, a better approximation can be found by using the Taylor series of the square root

As an iterative method, the order of convergence is equal to the number of terms used. With two terms, it is identical to the Babylonian method. With three terms, each iteration takes almost as many operations as the Bakhshali approximation, but converges more slowly. Therefore, this is not a particularly efficient way of calculation. To maximize the rate of convergence, choose is as small as possible.

Continued fraction expansion

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 irrationals, where a, b and c are integers), and in particular, square roots of integers, have periodic continued fractions. 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, we can express the square root of S

Proceeding this way, we get a generalized continued fraction

The first step to evaluating such a fraction 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, and for ,

rounded to three digits of precision.

The actual value of , good to almost 4 digits of precision, etc.

The following are examples of square roots, their simple continued fractions, and their first terms called convergents up to and including denominator 99: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 e.g., no fraction with a denominator less than or equal to 70 is as good an approximation to as 99/70.

Approximations that depend on the floating point representation

A number is represented in a floating point which is also called scientific notation is good to an order of magnitude. Next, recognise that some powers, , will be odd, thus for 3141.59 = 3.14159 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 so that the square root will be .

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 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 22 an even power so the paired bits of the mantissa are 01, while .625 = 0.101 in binary normalises to 1.01 2−1 an odd power so the adjustment is to 10.1 2−2 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 EXPONENT(x) 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

So for a 32-bit single precision floating point number in IEEE format (where notably, the power has a bias

For example, 1.0 is represented by a hexadecimal number . In a similar fashion you get 0.5 from 1.5 ().

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 is the exponent bias and

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

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 = , the maximum relative error is minimized to ±3.5%. For a = 0, the approximation is greater than or equal to the square root of val for any value of val.

If the approximation is to be used for an initial guess for Newton's method, then the reciprocal form shown in the following section is preferred.

The approximation can be further optimized by combining the addition of a, 1 << 29 and -1 << 22 into a single operation. This results in (val.i >> 1) + 0x1FBB4F2E for a = that minimizes the error correction and (val.i >> 1) + 1FC00000 for a = 0.

Reciprocal of the square root

A variant of the above routine is included below, which can be used to compute the reciprocal 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. In computer graphics it is a very efficient way to normalize a vector.

Some VLSI hardware implements inverse square root using a second degree polynomial estimation followed by a Goldschmidt iteration.

Negative or complex square

If S<

If S=a+bi where a and b are real and b≠

This can be verified by squaring the root.

is the 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

• Integer square root

• Alpha max plus beta min algorithm

• nth root algorithm

• Fast inverse square root

Notes

References

Bibliography

External links