What algorithm do computers use to compute the square root of a number?

Question

What algorithm do computers use to compute the square root of a number ?

EDIT

It seems there is a similar question here: Finding square root without division and initial guess

But I like the answers provided here more. Plus the person asking the similar question had the question formed in a personal way.

My question has an easier to find wording. The similar one was asked having an insight about the inner mechanics of such an algorithm and was being asked like if it was something that had not yet been solved-found.

I don't know if there is a way for a moderator to merge the two posts.

Answer 1

Have a look at the paper

https://www.lirmm.fr/arith18/papers/burgessn-divider.pdf

which describes the division / square root hardware implementation in some ARM processors.

Roughly the idea: For a division unit, you can look up the initial bits of divisor and dividend in a table to get two bits of result. For double precision with a 53 bit mantissa, you repeat that 26 or 27 times and you have the result. For single precision with 24 bit mantissa, you need to repeat only 12 or 13 times. All the other bits, calculating the exponent etc., can be done in parallel.

Now what does square root have to do with division? Simple: sqrt (x) = x / sqrt (x). You get the square root by dividing x by its square root. Now when you start this division, you don't know the square root yet. But as soon as you have say eight bits and you are very, very careful, you can produce two mantissa bits per cycle. And getting the first say eight bits is not hard, since unlike a division x / y, there is only one operand involved.

The paper describes this in much more detail, and importantly describes why this is superior for example to using Newton-Raphson or similar. Note: It is superior for typical floating-point numbers with 24 or 53 bit mantissa. Since that is what is used 99.9999% of the time it is important, but for different situations things are different.

Answer 2

Well, this is kind of a numerical method question and it is not related to approximation algorithm. Anyway I think you can search for Bisection method, Newton's method (or Newton-Raphason method), there are of course other ways but this is what I remember from a course in numerical methods.

Check this textbook:

Steven C. Chapra, and Ryamond P. Canale. Numerical Methods. 7th edition. McGraw Hill Education.

Answer 3

If you are not a hardware designer, then you likely use the Newton iteration method. Given an equation f(x) = 0, and an approximate solution $x_0$, Newton iteration calculates a (hopefully) better approximation as $x_{n+1} = x_n - f(x_n) / f'(x_n)$. If we change $x = a^{1/2}$ to $x^2 = a$ to $x^2 - a = 0$, then this gives the simple formula $x_{n+1}$ = $x_n - ({x_n}^2 - a) / 2x_n$ = $(x_n + a/x_n)/2$.

This works well with a good initial approximation, but each iteration round includes a division, and divisions are sloooow. We therefore use a different formula: Instead of solving $x = a^{1/2}$ we solve $x = a^{-1/2}$ and multiply the result by a, which gives $a^{1/2}$ as requested.

We rearrange $x = a^{-1/2}$ as $x^2 = 1/a$, then $1/x^2 = a$ and $1/x^2 - a = 0$. Now $f'(x) = -2/x^3$, so Newton-iteration gives $x_n - (1/{x_n}^2 - a) / (-2/{x_n}^3)$ or $x_n + (x_n - a{x_n}^3) / 2)$ or $1.5x_n - (a/2){x_n}^3))$. This can be calculated using multiplications only, which is much faster than division. On a modern processor which can calculate multiplications in parallel and has a fused multiply-add operation, you calculate $r = 1.5x_n$, $s = (a/2)x_n$, $t = {x_n}^2$, then $x_{n+1} = r + s \cdot t$.

To support this, x86 processors for example have a very fast hardware instruction calculating an approximation to $a^{-1/2}$ with a relative error less than 1 / 1024, implemented using a table.

Answer 4

You get a completely different answer if you ask for square roots of n-bit floating point numbers for LARGE n, say n ≥ 1000 or n ≥ 1,000,000.

For large n, we usually define the time it takes to multiply two n-bit numbers as M(n). That time depends on how clever our algorithm is; what we learned in school has M(n) = O (n^2), the Karatsuba algorithm has M(n) ≈ O (n^1.58), and FFT-based algorithms run in O (n log n). Division with n bit precision can be done in not much more time, about 2 M(n).

We can use the Newton-Raphson method, iterating $x_{k+1}$ = $x_k - ({x_k}^2 - a) / 2x_k$. However, we make one critical change: If we want a result with say 1,000,000 bits precision, we don't do the whole calculation with that precision. We calculate one approximation with just enough precision so that the final iteration step will give the required precision. For example, we might calculate the second-to-last value with 500,005 or 500,010 bits of precision and then do the last iteration with full precision. But calculating ${x_k}^2$, we know that the first 500,000 bits of the result should be identical with a, so we can save a lot of time not calculating those 500,000 bits when we already know what they are. And we don't need to calculate $({x_k}^2 - a) / 2x_k$ with 1,000,000 bit precision because the result is only a tiny correction to $x_k$ so just about 500,000 bits precision are enough. (Obviously you need to analyse error bounds very carefully).

When you calculate the second-to-last approximation, you also don't need 500,010 bits of precision throughout the whole calculation; it is enough to calculate the third-to-last approximation with say 250,015 bit precision and so on. So the total time is O(M(n)) + O(M(n/2)) + O(M(n/4)) + .... The cost of the last step dominates the rest as long as M(n) ≥ O(n), so a square root with n bit precision can be calculated in time c * M(n) for c probably around 6 to 10.

If you do many such calculations, it would be worthwhile checking whether calculating an approximation to $a^{-1/2}$ would give faster results.

Answer 5

Thanks to @YOUSEFY in the textbook of

Steven C. Chapra, and Ryamond P. Canale. Numerical Methods. 7th edition. McGraw Hill Education.

I found an iterative method called "divide and average" where a is the number in the root.

x = (x+a/x)/2

But what I still don't know is if this is the preferred way to compute it (in terms of optimization).