From Wikipedia, the free encyclopedia  View original article

In numerical analysis, a branch of mathematics, there are several square root algorithms or methods for calculating the principal square root of a nonnegative real number. For the square roots of a negative or complex number, see below.
Finding is the same as solving the equation . Therefore, any general numerical rootfinding algorithm can be used. Newton's method, for example, reduces in this case to the socalled Babylonian method:
Generally, these methods yield approximate results. To get a higher precision for the root, a higher precision for the square is required and a larger number of steps must be calculated.
Many square root algorithms require an initial seed value. If the initial seed value is far away from the actual square root, the algorithm will be slowed down. It is therefore useful to have a rough estimate, which may be very inaccurate but easy to calculate. With expressed in scientific notation as where and n is an integer, the square root can be estimated as
The factors two and six are used because they approximate the geometric means of the lowest and highest possible values with the given number of digits: and .
For , the estimate is .
When working in the binary numeral system (as computers do internally), by expressing as where , the square root can be estimated as , since the geometric mean of the lowest and highest possible values is .
For the binary approximation gives
These approximations are useful to find better seeds for iterative algorithms, which results in faster convergence.
Perhaps the first algorithm used for approximating is known as the "Babylonian method", named after the Babylonians,^{[1]} or "Hero's method", named after the firstcentury Greek mathematician Hero of Alexandria who gave the first explicit description of the method.^{[2]} It can be derived from (but predates by 16 centuries) Newton's method. The basic idea is that if x is an overestimate to the square root of a nonnegative real number S then will be an underestimate and 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 roots, thus assuring convergence).
More precisely, if is our initial guess of and is the error in our estimate such that , then we can expand the binomial and solve for
Therefore, we can compensate for the error and update our old estimate as
This becomes our next best guess. The process of updating is iterated until desired accuracy is obtained. This is a quadratically convergent algorithm, which means that the number of correct digits of the approximation roughly doubles with each iteration. It proceeds as follows:
It can also be represented as:
This algorithm works equally well in the padic numbers, but cannot be used to identify real square roots with padic 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 2adics.
To calculate , where S = 125348, to 6 significant figures, use the rough estimation method above to get
Therefore,
Let the relative error in x_{n} be defined by
and thus
Then it can be shown that
and thus that
and consequently that convergence is assured provided that x_{0} and S are both positive.
If using the rough estimate above with the Babylonian method, then the worst cases are^{[clarification needed]}:
Thus in any case,
Remember that rounding errors will slow the convergence. It is recommended to keep at least one extra digit beyond the desired accuracy of the x_{n} being calculated to minimize round off error.
This is a method to find each digit of the square root in a sequence. It is slower than the Babylonian method (if you have a calculator that can divide in one operation), but it has several advantages:
Napier's bones include an aid for the execution of this algorithm. The shifting nthroot algorithm is a generalization of this method.
Suppose we are to find the square root of N by expressing it as a sum of n positive numbers such that
By repeatedly applying the basic identity
the righthandside term can be expanded as
This expression allows us to find the values of s by making suitable guesses in a sequential manner. Each new guess should satisfy the recursion
where and Here at each recursion, the value of is to be found such that When the exact square root has been found; if not, then the sum of s gives a suitable approximation of the square root, with being the approximation error.
Suppose we desire to find three digit square root for a five digit number such that . Here the 1, 10, and 100 are merely the place holders for a three digit decimal number. In each recursion the value of is searched from the set of decimal digits
1. Find the first digit such that comes as close to as possible. That is,
2. Next, find the second digit such that comes as close to as possible. That is,
3. Finally, find the last digit such that comes as close to as possible. That is,
The if is an exact square. If not, then we can increase the precision by adding more terms in the decimal expansion.
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 digitbydigit square root in the binary number system is quite efficient since the value of is searched from a smaller set of binary digits {0,1}. This makes the computation faster at each stage since the value of is either zero for or for . Also the decision process is easier since all we need to do is check for if this number is greater than or less than to use it at the mth stage of calculation. Also, the fact that multiplication by 2 is done by left bitshifts helps in the computation.
Write the original number in decimal form. The numbers are written similar to the long division algorithm, and, as in long division, the root will be written on the line above. Now separate the digits into pairs, starting from the decimal point and going both left and right. The decimal point of the root will be above the decimal point of the square. One digit of the root will appear above each pair of digits of the square.
Beginning with the leftmost pair of digits, do the following procedure for each pair:
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 is 12.34
Find the square root of 2.
1. 4 1 4 2 / \/ 02.00 00 00 00 02 1*1 <= 2 < 2*2 x = 1 01 y = x*x = 1*1 = 1 01 00 24*4 <= 100 < 25*5 x = 4 00 96 y = (20+x)*x = 24*4 = 96 04 00 281*1 <= 400 < 282*2 x = 1 02 81 y = (280+x)*x = 281*1 = 281 01 19 00 2824*4 <= 11900 < 2825*5 x = 4 01 12 96 y = (2820+x)*x = 2824*4 = 11296 06 04 00 28282*2 <= 60400 < 28283*3 x = 2 The desired precision is achieved: The square root of 2 is about 1.4142
Inherent to digitbydigit algorithms is a search and test step: find a digit, , when added to the right of a current solution , such that , where is the value for which a root is desired. Expanding: . The current value of —or, usually, the remainder—can be incrementally updated efficiently when working in binary, as the value of will be a single bit, and the operations needed to compute and can be replaced with faster bit shift operations.
Here we obtain the square root of 81, which when converted into binary gives 1010001. The numbers in the left column gives the option between that number or zero to be used for subtraction at that stage of computation. The final answer is 1001, which in decimal is 9.
1 0 0 1  √ 1010001 1 1 1  101 01 0  1001 100 0  10001 10001 10001  0
This gives rise to simple computer implementations:^{[3]}
short isqrt(short num) { short res = 0; short bit = 1 << 14; // The secondtotop bit is set: 1L<<30 for long // "bit" starts at the highest power of four <= the argument. while (bit > num) bit >>= 2; while (bit != 0) { if (num >= res + bit) { num = res + bit; res = (res >> 1) + bit; } else res >>= 1; bit >>= 2; } return res; }
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.^{[4]}
Pocket calculators 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 () and exponentials ():
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 tables or slide rules.
This method for finding an approximation to a square root was described in an ancient Indian mathematical manuscript called the Bakhshali manuscript. It is equivalent to two iterations of the Babylonian method beginning with N. The original presentation goes as follows: To calculate , let N^{2} be the nearest perfect square to S. Then, calculate:
This can be also written as:
Find
The Vedic duplex method from the book 'Vedic mathematics' is a variant of the digit by digit method for calculating the square root.^{[5]} The duplex is the square of the central digit plus double the crossproduct of digits equidistant from the center. The duplex is computed from the quotient digits (square root digits) computed thus far, but after the initial digits. The duplex is subtracted from the dividend digit prior to the second subtraction for the product of the quotient digit times the divisor digit. For perfect squares the duplex and the dividend will get smaller and reach zero after a few steps. For nonperfect squares the decimal value of the square root can be calculated to any precision desired. However, as the decimal places proliferate, the duplex adjustment gets larger and longer to calculate. The duplex method follows the Vedic ideal for an algorithm, oneline, mental calculation. It is flexible in choosing the first digit group and the divisor. Small divisors are to be avoided by starting with a larger initial group.
In short, to calculate the duplex of a number, double the product of each pair of equidistant digits plus the square of the center digit (of the digits to the right of the colon).
Number => Calculation = Duplex 574 ==> 2(5·4) + 7^{2} = 89 406,739 ==> 2(4·9)+ 2(0·3)+ 2(6·7) = 72+0+84 = 156 123,456 ==> 2(1·6)+ 2(2·5)+ 2(3·4) = 12 +20 +24 = 56 88,900,777 ==> 2(8·7)+2(8·7)+2(9·7)+2(0·0) = 112+112+126+0 = 350 48329,03711 ==> 2(4·1)+2(8·1)+2(3·7)+2(2·3)+2(9·0)= 8+16+42+12+0 = 78
In a square root calculation the quotient digit set increases incrementally for each step.
Number => Calculation = Duplex: 1 ==> 1^{2} = 1 14 ==>2(1·4) = 8 142 ==> 2(1·2) + 4^{2} = 4 + 16 = 20 14,21 ==> 2(1·1) + 2(4·2) = 2 + 16 = 18 14213 ==> 6+8+4 = 18 142,135 ==> 10+24+4 = 38 1421356 ==> 12+40+12+1 = 65 1421,3562 ==> 4+48+20+6 = 78 142,135,623 ==> 6+16+24+10+9 = 65 142,1356,237 ==> 14+24+8+12+30 = 88 142,13562,373 ==> 6+56+12+4+36+25 = 139
Consider the perfect square 2809 = 53^{2}. Use the duplex method to find the square root of 2,809.
Find the square root of 2809. Set down the number in groups of two digits. The number of groups gives the number of whole digits in the root. Put a colon after the first group, 28, to separate it. From the first group, 28, obtain the divisor, 10, since 28>25=5^{2} and by doubling this first root, 2x5=10. Gross dividend: 28: 0 9. Using mental math: Divisor: 10) 3 0 Square: 10) 28: _{3}0 9 Duplex, Deduction: 25: xx 09 Square root: 5: 3. 0 Dividend: 30 00 Remainder: 3: 00 00 Square Root, Quotient: 5: 3. 0
Find the square root of 2,080,180,881. Solution by the duplex method: this tendigit square has five digitpairs, so it will have a fivedigit square root. The first digitpair is 20. Put the colon to the right. The nearest square below 20 is 16, whose root is 4, the first root digit. So, use 2·4=8 for the divisor. Now proceed with the duplex division, one digit column at a time. Prefix the remainder to the next dividend digit.
divisor; gross dividend: 8) 20: 8 0 1 8 0 8 8 1 read the dividend diagonally up: 4 8 7 11 10 10 0 8 minus the duplex: 16: xx 25 60 36 90 108 00 81 actual dividend: : 48 55 11 82 10 00 08 00 minus the product: : 40 48 00 72 00 00 0 00 remainder: 4: 8 7 11 10 10 0 8 00 quotient: 4: 5, 6 0 9. 0 0 0 0
Duplex calculations: Quotientdigits ==> Duplex deduction. 5 ==> 5^{2}= 25 5 and 6 ==> 2(5·6) = 60 5,6,0 ==> 2(5·0)+6^{2} = 36 5,6,0,9 ==> 2(5·9)+2(6·0) = 90 5,6,0,9,0 ==> 2(5·0)+2(6·9)+ 0 = 108 5,6,0,9,0,0 ==> 2(5·0)+2(6·0)+2(0·9) = 0 5,6,0,9,0,0,0 ==> 2(5·0)+2(6·0)+2(0·0)+9^{2} = 81
Hence the square root of 2,080,180,881 is exactly 45,609.
Find the square root of two to ten places. Take 20,000 as the beginning group, using three digitpairs at the start. The perfect square just below 20,000 is 141, since 141^{2} = 19881 < 20,000. So, the first root digits are 141 and the divisor doubled, 2 x 141 = 282. With a larger divisor the duplex will be relatively small. Hence, the multiple of the divisor can be picked without confusion.
Dividend: 2.0000 : 0 0 0 0 0 0 0 0 Diagonal;Divisor: (282) : 1190 620 400 1020 1620 1820 750 1120 Minus duplex: : xxxx 16 16 12 28 53 74 59 Actual dividend: 20000 : 1190 604 384 1008 1592 1767 676 1061 Minus product: 19881 : 1128 564 282 846 1410 1692 564 846 Remainder: 119 : 62 40 102 162 182 75 112 215 Root quotient: 1.41 : 4 2 1 3 5 6 2 3
Ten multiples of 282: 282; 564; 846; 1128; 1410; 1692; 1974; 2256; 2538; 2820.
This method is applicable for finding the square root of and converges best for . 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 by an integer power of 4, and therefore by the corresponding power of 2, by changing the exponent or by shifting, respectively. Therefore, 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 is
while the iterative steps read
Then, (while ).
Note that the convergence of , and therefore also of , is quadratic.
The proof of the method is rather easy. First, rewrite the iterative definition of as
Then it is straightforward to prove by induction that
and therefore the convergence of to the desired result is ensured by the convergence of to 0, which in turn follows from .
This method was developed around 1950 by M. V. Wilkes, D. J. Wheeler and S. Gill^{[6]} for use on EDSAC, one of the first electronic computers.^{[7]} The method was later generalized, allowing the computation of nonsquare roots.^{[8]}
The following are iterative methods for finding the reciprocal square root of S which is . Once it has been found, find by simple multiplication: . 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.
Some computers use Goldschmidt's algorithm to simultaneously calculate and . Goldschmidt's algorithm finds faster than NewtonRaphson iteration on a computer with a fused multiply–add instruction and either a pipelined floating point unit or two independent floatingpoint units. Two ways of writing Goldschmidt's algorithm are:^{[9]}
Each iteration:
until is sufficiently close to 1, or a fixed number of iterations.
which causes
Goldschmidt's equation can be rewritten as:
Each iteration: (All 3 operations in this loop are in the form of a fused multiply–add.)
until is sufficiently close to 0, or a fixed number of iterations.
which causes
If N is an approximation to , a better approximation can be found by using the Taylor series of the square root function:
As an iterative method, the order of convergence is equal to the number of terms used. With 2 terms, it is identical to the Babylonian method; With 3 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 N so that is as small as possible.
Other methods are less efficient than the ones presented above.
A completely different method for computing the square root is based on the CORDIC algorithm, which uses only very simple operations (addition, subtraction, bitshift and table lookup, but no multiplication). The square root of S may be obtained as the output using the hyperbolic coordinate system in vectoring mode, with the following initialization:^{[10]}
Quadratic irrationals (numbers of the form , 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. The following iterative algorithm can be used for this purpose (S is any natural number that is not a perfect square):
Notice that m_{n}, d_{n}, and a_{n} are always integers. The algorithm terminates when this triplet is the same as one encountered before. The algorithm can also terminate on a_{i} when a_{i} = 2 a_{0},^{[11]} which is easier to implement.
The expansion will repeat from then on. The sequence [a_{0}; a_{1}, a_{2}, a_{3}, …] is the continued fraction expansion:
Begin with m_{0} = 0; d_{0} = 1; and a_{0} = 10 (10^{2} = 100 and 11^{2} = 121 > 114 so 10 chosen).
So, m_{1} = 10; d_{1} = 14; and a_{1} = 1.
Next, m_{2} = 4; d_{2} = 7; and a_{2} = 2.
Now, loop back to the second equation above.
Consequently, the simple continued fraction for the square root of 114 is
Its actual value is approximately 10.67707 82520 31311 21....
A more rapid method is to evaluate its generalized continued fraction. From the formula derived there:
and the fact that 114 is 2/3 of the way between 10^{2}=100 and 11^{2}=121 results in
which is simply the aforementioned [10;1,2, 10,2,1, 20,1,2, 10,2,1, 20,1,2, ...] evaluated at every third term. Combining pairs of fractions produces
which is now [10;1,2, 10,2,1,20,1,2, 10,2,1,20,1,2, ...] evaluated at the third term and every six terms thereafter.
Pell's equation (also known as Brahmagupta equation since he was the first to give a solution to this particular equation) and its variants yield a method for efficiently finding continued fraction convergents of square roots of integers. However, it can be complicated to execute, and usually not every convergent is generated. The ideas behind the method are as follows:
The method is as follows:
On computers, a very rapid Newton'smethodbased approximation to the square root can be obtained for floating point numbers when computers use an IEEE (or sufficiently similar) representation.
This technique is based on the fact that the IEEE floating point format approximates base2 logarithm. For example, you can get the approximate logarithm of 32bit single precision floating point number by translating its binary representation as an integer, scaling it by , and removing a bias of 127, i.e.
For example, 1.0 is represented by a hexadecimal number 0x3F800000, which would represent if taken as an integer. Using the formula above you get , as expected from . In a similar fashion you get 0.5 from 1.5 (0x3FC00000).
To get the square root, divide the logarithm by 2 and convert the value back. The following program demonstrates the idea. Note that 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 is the number of explicitly stored bits in the mantissa and then show that
/* Assumes that float is in the IEEE 754 single precision floating point format * and that int is 32 bits. */ float sqrt_approx(float z) { int val_int = *(int*)&z; /* Same bits, but as an int */ /* * To justify the following code, prove that * * ((((val_int / 2^m)  b) / 2) + b) * 2^m = ((val_int  2^m) / 2) + ((b + 1) / 2) * 2^m) * * where * * b = exponent bias * m = number of mantissa bits * * . */ val_int = 1 << 23; /* Subtract 2^m. */ val_int >>= 1; /* Divide by 2. */ val_int += 1 << 29; /* Add ((b + 1) / 2) * 2^m. */ return *(float*)&val_int; /* Interpret again as float */ }
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
val_int = (1 << 29) + (val_int >> 1)  (1 << 22) + a;
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 = 0x4C000, the errors are between about 3.5% and 3.5%.
If the approximation is to be used for an initial guess for Newton's method to the equation , then the reciprocal form shown in the following section is preferred.
A variant of the above routine is included below, which can be used to compute the reciprocal of the square root, i.e., instead, was written by Greg Walsh, and implemented into SGI Indigo by Gary Tarolli.^{[12]}^{[13]} The integershift 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.^{[14]} In computer graphics it is a very efficient way to normalize a vector.
float invSqrt(float x) { float xhalf = 0.5f*x; union { float x; int i; } u; u.x = x; u.i = 0x5f3759df  (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; }
Some VLSI hardware implements inverse square root using a second degree polynomial estimation followed by a Goldschmidt iteration.^{[15]}
If S < 0, then its principal square root is
If S = a+bi where a and b are real and b ≠ 0, then its principal square root is
This can be verified by squaring the root.^{[16]}^{[17]} Here
is the modulus of S. The principal square root of a complex number is defined to be the root with the nonnegative real part.