\sqrt{S}
S
Most square root computation methods are iterative: after choosing a suitable initial estimate of
\sqrt{S}
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).
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 "digits" after the 1, but it is not known exactly how. They knew how to approximate a hypotenuse using(giving for example
41 | + | |
60 |
15 | |
3600 |
40 | |
60 |
10 | |
60 |
\sqrt2.
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.
Many iterative square root algorithms require an initial seed value. The seed must be a non-zero positive number; it should be between 1 and
S
x0=1
S
\tfrac12\vertlog2S\vert
In general, an estimate is pursuant to an arbitrary interval known to contain the root (such as
[x0,S/x0]
f(x)=\sqrt{x}
f(x)
Typically the number
S
a x 102n
1\leqa<100
\sqrta x 10n
1\leq\sqrta<10
Scalar methods divide the range into intervals, and the estimate in each interval is represented by a single scalar number. If the range is considered as a single interval, the arithmetic mean (5.5) or geometric mean (
\sqrt{10} ≈ 3.16
10n
For two intervals, divided geometrically, the square root
\sqrt{S}=\sqrt{a} x 10n
This estimate has maximum absolute error of
4 ⋅ 10n
For example, for
S=125348
12.5348 x 104
\sqrt{S} ≈ 6 ⋅ 102=600
\sqrt{125348}=354.0
A better estimate, and the standard method used, is a linear approximation to the function
y=x2
S
[1,100]
A least-squares regression line minimizes the average difference between the estimate and the value of the function. Its equation is
y=8.7x-10
x=0.115y+1.15
That is the best estimate on average that can be achieved with a single piece linear approximation of the function y=x2 in the interval
[1,100]
To divide by 10, subtract one from the exponent of
a
[1,100]
A much better estimate can be obtained by a piece-wise linear approximation: multiple line segments, each approximating some subarc of the original. The more line segments used, the better the approximation. The most common way is to use tangent lines; the critical choices are how to divide the arc and where to place the tangent points. An efficacious way to divide the arc from y = 1 to y = 100 is geometrically: for two intervals, the bounds of the intervals are the square root of the bounds of the original interval, 1×100, i.e. [1,{{radic|100|2}}] and [{{radic|100|2}},100]. For three intervals, the bounds are the cube roots of 100: [1,{{radic|100|3}}], [{{radic|100|3}},({{radic|100|3}})<sup>2</sup>], and [({{radic|100|3}})<sup>2</sup>,100], etc. For two intervals, = 10, a very convenient number. Tangent lines are easy to derive, and are located at x = and x = . Their equations are:
y=3.56x-3.16
y=11.2x-31.6
x=0.28y+0.89
x=.089y+2.8
S=a ⋅ 102n
The maximum absolute errors occur at the high points of the intervals, at a=10 and 100, and are 0.54 and 1.7 respectively. The maximum relative errors are at the endpoints of the intervals, at a=1, 10 and 100, and are 17% in both cases. 17% or 0.17 is larger than 1/10, so the method yields less than a decimal digit of accuracy.
In some cases, hyperbolic estimates may be efficacious, because a hyperbola is also a convex curve and may lie along an arc of Y = x2 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
[1,100]
S=a ⋅ 102n
The floating division need be accurate to only one decimal digit, because the estimate overall is only that accurate, and can be done mentally. A hyperbolic estimate is better on average than scalar or linear estimates. It has maximum absolute error of 1.58 at 100 and maximum relative error of 16.0% at 10. For the worst case at a=10, the estimate is 3.67. If one starts with 10 and applies Newton-Raphson iterations straight away, two iterations will be required, yielding 3.66, before the accuracy of the hyperbolic estimate is exceeded. For a more typical case like 75, the hyperbolic estimate is 8.00, and 5 Newton-Raphson iterations starting at 75 would be required to obtain a more accurate result.
A method analogous to piece-wise linear approximation but using only arithmetic instead of algebraic equations, uses the multiplication tables in reverse: the square root of a number between 1 and 100 is between 1 and 10, so if we know 25 is a perfect square (5 × 5), and 36 is a perfect square (6 × 6), then the square root of a number greater than or equal to 25 but less than 36, begins with a 5. Similarly for numbers between other squares. This method will yield a correct first digit, but it is not accurate to one digit: the first digit of the square root of 35 for example, is 5, but the square root of 35 is almost 6.
A better way is to the divide the range into intervals 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.[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:
nearest square | k=\sqrt{a} | ||
---|---|---|---|
1 to 2.5 | 1 (= 12) | 1 | |
2.5 to 6.5 | 4 (= 22) | 2 | |
6.5 to 12.5 | 9 (= 32) | 3 | |
12.5 to 20.5 | 16 (= 42) | 4 | |
20.5 to 30.5 | 25 (= 52) | 5 | |
30.5 to 42.5 | 36 (= 62) | 6 | |
42.5 to 56.5 | 49 (= 72) | 7 | |
56.5 to 72.5 | 64 (= 82) | 8 | |
72.5 to 90.5 | 81 (= 92) | 9 | |
90.5 to 100 | 100 (= 102) | 10 |
The final operation is to multiply the estimate by the power of ten divided by 2, so for
S=a ⋅ 102n
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
k2\lea<(k+1)2
\sqrt{a}
The final operation, as above, is to multiply the result by the power of ten divided by 2;
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.
Example: find the square root of 75., so is 75 and is 0. From the multiplication tables, the square root of the mantissa must be 8 point something because 8 × 8 is 64, but 9 × 9 is 81, too big, so is 8; something is the decimal representation of . The fraction is 75 - k2 = 11, the numerator, and 81 - k2 = 17, the denominator. 11/17 is a little less than 12/18, which is 2/3s or .67, so guess .66 (it's ok to guess here, the error is very small). So the estimate is . to three significant digits is 8.66, so the estimate is good to 3 significant digits. Not all such estimates using this method will be so accurate, but they will be close.
When working in the binary numeral system (as computers do internally), by expressing
S
a x 22n
0.12\leqa<102
\sqrt{S}=\sqrt{a} x 2n
which is the least-squares regression line to 3 significant digit coefficients.
\sqrt{a}
a
a
\sqrt{S} ≈ (0.5+0.5 ⋅ a) ⋅ 2n
For
S=125348=1 1110 1001 1010 01002=
16 | |
1.1110 1001 1010 0100 | |
2 x 2 |
\sqrt{S} ≈ (0.5+0.5 ⋅ a) ⋅ 28=1.0111 0100 1101 00102 ⋅ 1 0000 00002=1.456 ⋅ 256=372.8
\sqrt{125348}=354.0
An estimate for
a
a
The first explicit algorithm for approximating
\sqrt{S~}
Given a positive real number
S
l( x0, x1, x2, x3, \ldots r)
\limnxn=\sqrt{S~}~.
This is equivalent to using Newton's method to solve
x2-S=0
xn
The basic idea is that if
x
S
\tfrac{ S }{x}
More precisely, if
x
\sqrt{S~}
\varepsilon
S=\left(x+\varepsilon\right)2 ,
\varepsilon=
S-x2 | |
2x+\varepsilon |
≈
S-x2 | |
2x |
,
\varepsilon\llx~
This algorithm works equally well in the -adic numbers, 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.
To calculate
\sqrt{S~}
S=125348 ,
Therefore
\sqrt{125348~} ≈ 354.045
Suppose that
x0>0~~and~~S>0~.
n:xn>0~.
xn
Then it can be shown that
And thus thatand consequently that convergence is assured, and quadratic.
If using the rough estimate above with the Babylonian method, then the least accurate cases in ascending order are as follows:
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
xn
This method for finding an approximation to a square root was described in an ancient Indian manuscript, called the Bakhshali manuscript. It is equivalent to two iterations of the Babylonian method beginning with x0. Thus, the algorithm is quartically convergent, which means that the number of correct digits of the approximation roughly quadruples with each iteration. The original presentation, using modern notation, is as follows: To calculate
\sqrt{S}
2 | |
x | |
0 |
S
This can be used to construct a rational approximation to the square root by beginning with an integer. If
x0=N
N2
S
d=S-N2
The Bakhshali method can be generalized to the computation of an arbitrary root, including fractional roots.
Using the same example as given with the Babylonian method, let
S=125348.
Likewise the second iteration gives
This is a method to find each digit of the square root in a sequence. This method is based on the binomial theorem and basically an inverse algorithm solving
(x+y)2=x2+2xy+y2
Disadvantages are:
Napier's bones include an aid for the execution of this algorithm. The shifting th root algorithm is a generalization of this method.
First, consider the case of finding the square root of a number, that is the square of a two-digit number, where is the tens digit and is the units digit. Specifically:
Now using the digit-by-digit algorithm, we first determine the value of . is the largest digit such that is less than or equal to from which we removed the two rightmost digits.
In the next iteration, we pair the digits, multiply by 2, and place it in the tenth's place while we try to figure out what the value of is.
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 positive numbers such that
By repeatedly applying the basic identity the right-hand-side term can be expanded as
This expression allows us to find the square root by sequentially guessing the values of
ai
a1,\ldots,am-1
Ym=\left[2Pm-1+am\right]am,
am
Xm\geq0
1\leqm\leqn,
X0=N.
Xn=0,
ai
Xn
For example, in the decimal number system we have where
10n-i
ai\in\{0,1,2,\ldots,9\}
Pm-1
Ym
Here since the place value of
Ym
Xm-1
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
ai
Ym
Ym=0
am=0
Ym=2Pm-1+1
am=1
am
am
Ym\leqXm-1
am=1.
am=1
am=0.
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:
x(20p+x)\lec
x
Find the square root of 152.2756.
1 2. 3 4 / \/ 01 52.27 56 01 1*1 <= 1 < 2*2 x=1 01 y = x*x = 1*1 = 1 00 52 22*2 <= 52 < 23*3 x=2 00 44 y = (20+x)*x = 22*2 = 44 08 27 243*3 <= 827 < 244*4 x=3 07 29 y = (240+x)*x = 243*3 = 729 98 56 2464*4 <= 9856 < 2465*5 x=4 98 56 y = (2460+x)*x = 2464*4 = 9856 00 00 Algorithm terminates: Answer=12.34
This section uses the formalism from the digit-by-digit calculation section above, with the slight variation that we let
N2=(an+...b+
2 | |
a | |
0) |
am=2m
am=0
2m
2n
20
Pm=an+an-1+\ldots+am
ai
am
2m
0
Pm=Pm+1+2m
2 | |
P | |
m |
\leqN2
2m
am=2m
am=0
Pm=Pm+1
Pm
Xm=N2-
2 | |
P | |
m |
Xm=Xm+1-Ym
Ym=
2 | |
P | |
m |
-
2 | |
P | |
m+1 |
=2Pm+1am+
2 | |
a | |
m |
an=Pn=2n
n
(2n)2=4n\leqN2
As an extra optimization, we store
Pm+12m+1
(2m)2
Ym
am
cm
dm
cm
dm
Note that: which is the final result returned in the function below.
An implementation of this algorithm in C:
Faster algorithms, in binary and decimal or any other base, can be realized by using lookup tables—in effect trading more storage space for reduced run time.
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 (
lnxn=nlnx
eln=x
This method is applicable for finding the square root of
0<S<3
S ≈ 1
S
\sqrt{S}
S
The initialization step of this method iswhile the iterative steps readThen,
an\to\sqrt{S}
cn\to0
The convergence of
cn
an
The proof of the method is rather easy. First, rewrite the iterative definition of
cn
an
\sqrt{S}
cn
-1<c0<2
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.
The following are iterative methods for finding the reciprocal square root of S which is
1/\sqrt{S}
\sqrt{S}
\sqrt{S}=S ⋅ (1/\sqrt{S})
(1/x2)-S=0
Goldschmidt's algorithm is an extension of Goldschmidt division, named after Robert Elliot Goldschmidt,[5] [6] which can be used to calculate square roots. Some computers use Goldschmidt's algorithm to simultaneously calculate
\sqrt{S}
1/\sqrt{S}
\sqrt{S}
The first way of writing Goldschmidt's algorithm begins
b0=S
Y0 ≈ 1/\sqrt{S}
y0=Y0
x0=Sy0
bi
xn
yn
xn=Syn
A second form, using fused multiply-add operations, begins
y0 ≈ 1/\sqrt{S}
x0=Sy0
h0=\tfrac12y0
ri
If N is an approximation to
\sqrt{S}
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 so that
|d| | |
N2 |
See also: Solving quadratic equations with continued fractions. The continued fraction representation of a real number can be used instead of its decimal or binary expansion and this representation has the property that the square root of any rational number (which is not already a perfect square) has a periodic, repeating expansion, similar to how rational numbers have repeating expansions in the decimal notation system.
Quadratic irrationals (numbers of the form
a+\sqrt{b | |
S=a2+r.
S-a2=(\sqrt{S}+a)(\sqrt{S}-a)=r
By applying this expression for
\sqrt{S}
The numerator/denominator expansion for continued fractions (see left) is cumbersome to write as well as to embed in text formatting systems. So mathematicians have devised several alternative notations, like[7]
When
r=1
For, the value of
a
Proceeding this way, we get a generalized continued fraction for the square root as
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,
r
a
Step 2 is to reduce the continued fraction from the bottom up, one denominator at a time, to yield a rational fraction whose numerator and denominator are integers. The reduction proceeds thus (taking the first three denominators):
Finally (step 3), divide the numerator by the denominator of the rational fraction to obtain the approximate value of the root: rounded to three digits of precision.
The actual value of is 1.41 to three significant digits. The relative error is 0.17%, so the rational fraction is good to almost three digits of precision. Taking more denominators gives successively better approximations: four denominators yields the fraction
41 | |
29 |
=1.4137
The following are examples of square roots, their simple continued fractions, and their first terms called convergents up to and including denominator 99:
~decimal | continued fraction | convergents | ||||||||||||||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
1.41421 | [1;\overline{2}] |
,
,
,
,
| ||||||||||||||||||||||
1.73205 | [1;\overline{1,2}] |
,
,
,
,
,
,
| ||||||||||||||||||||||
2.23607 | [2;\overline{4}] |
,
,
| ||||||||||||||||||||||
2.44949 | [2;\overline{2,4}] |
,
,
,
| ||||||||||||||||||||||
3.16228 | [3;\overline{6}] |
,
| ||||||||||||||||||||||
\sqrt{\pi} | 1.77245 | [1;1,3,2,1,1,6...] |
,
,
,
,
| |||||||||||||||||||||
\sqrt{e} | 1.64872 | [1;1,1,1,5,1,1...] |
,
,
,
,
,
| |||||||||||||||||||||
\sqrt{\phi} | 1.27202 | [1;3,1,2,11,3,7...] |
,
,
|
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.
A number is represented in a floating point format as
m x bp
\sqrt{m} x bp/2
bp/2
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 (or sufficiently similar) representation, and a very rapid approximation to the square root can be obtained for starting Newton's method. The technique that follows is based on the fact that the floating point format (in base two) approximates the base-2 logarithm. That is
log2(m x 2p)=p+log2(m)
So for a 32-bit single precision floating point number in IEEE format (where notably, the power has a bias of 127 added for the represented form) you can get the approximate logarithm by interpreting its binary representation as a 32-bit integer, scaling it by
2-23
For example, 1.0 is represented by a hexadecimal number 0x3F800000, which would represent
1065353216=127 ⋅ 223
1065353216 ⋅ 2-23-127=0
log2(1.0)
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
b
n
float sqrt_approx(float z)
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 = -0x4B0D2, the maximum relative error is minimized to ±3.5%.
If the approximation is to be used for an initial guess for Newton's method to the equation
(1/x2)-S=0
See main article: Fast inverse square root. A variant of the above routine is included below, which can be used to compute the reciprocal of the square root, i.e.,
x-1/2
Some VLSI hardware implements inverse square root using a second degree polynomial estimation followed by a Goldschmidt iteration.
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.Here
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.
\sqrt{\sqrt{1} ⋅ \sqrt{10}}=\sqrt[4]{10} ≈ 1.78
\sqrt{\sqrt{10} ⋅ \sqrt{100}}=\sqrt[4]{1000} ≈ 5.62