# Square Root Algorithm

I would like an efficient algorithm for square root of a positive integer. Is there a reference that compares various square root algorithms? Thanks.

-

## migrated from mathoverflow.netJul 8 '13 at 19:39

This question came from our site for professional mathematicians.

Not sure this really qualifies as a research level question, you might be better off posting on MathStackexchage. But a very fast and easy algorithm to compute square root of $A$ to $N$ decimal places for some reasonable $N$ is Newton's algorithm to find a root of $X^2 - A$. So pick an initial $x_0$, say $x_0=A/2$, and then iteratively compute $x_{i+1} = \dfrac{x_i}{2}-\dfrac{A}{2x_i}$. – Joe Silverman Jul 8 '13 at 19:02
The documentation and source code to GMP could be used as a reference. Also, do note that the size of the integer is important; different methods are better for differently sized integers. Any other special properties the integer might have (such as being a perfect square, or near a power of 2) are also very relevant. – Hurkyl Jul 8 '13 at 19:42
Do you want the integer square root? That is, the integer $i$ such that $i^2 \le n \lt (i+1)^2$? Or do you want an approximate rational number $r$ with $r^2\approx n$? And do you want an algorithm that is suitable for pencil-and-paper, for an electronic digital computer, or for something else? – MJD Jul 8 '13 at 19:42
– lhf Jan 4 at 12:27

Method 1

Newton's Method converges quadratically. $a \ge 0, x_0 \ne \sqrt{a}$

$$x_{k+1} = \dfrac{1}{2}\left(x_{k} +\dfrac{a}{x_{k}}\right)$$

Method 2

$n \ge 0, x_n \gt \sqrt{a}$ for all $n \gt 0$.

$$x^2_{k+1}-a = \left[\dfrac{x^2_{k} - a}{2x_n}\right]^2$$

Method 3

Third order method, $n \ge 0$

$$x_{n+1} = \dfrac{x_n(x_n^2 + 3a)}{2x^2_n+a}$$

Method 4

-
Your method 1 is missing a factor of $1/2$, I think. – Daryl Jul 8 '13 at 22:13
@Daryl: Thank you for catching that, I accidentally deleted it before posting when I was cleaning up! Thanks! Regards – Amzoti Jul 8 '13 at 22:38
Nice details & overview (& links, to boot!) (+) – amWhy Jul 9 '13 at 0:06
@amWhy: Thank you. Surprising how many interesting algorithms there are for that! ;-) – Amzoti Jul 9 '13 at 1:19

The Newton's Method will not work if you are limited to integer arithmetic, which this problem seems to presuppose.

You know for sure that $s = \sqrt{N} \in [N/2, N]$.

Start with $x = N/2$ (integer division by 2) and check:

• if $x^2 = N$, you are done;
• if $x^2 < N$, you know that $s \in [x,N]$;
• if $x^2>N$ then you know $s \in [N/2,x]$.

Set $s$ to be the midpoint of the resulting interval and repeat the previous step.

The idea is to split the interval in half each iteration. The running time is worst-case $\log(N)$.

-
If $N \ge 9$, it is not true that $\sqrt N \in [N/2,N]$ – Ross Millikan Jul 8 '13 at 20:06
$\log(N)$ arithmetic operations is not $\log(N)$ time.... – Hurkyl Jul 8 '13 at 20:52

"efficient" rather depends on what you constraints are. For instance, if you have enough memory to store some floats, but little cpu time, then you can store a look up table for all integers until some power of 10. So say you had to find $\sqrt(1632397825)$. You could write: $$\sqrt{1732397825} =\sqrt{ 17*10^8 + 32397825}\sim \sqrt{17}\cdot 10^4 + \epsilon$$

To calculate $\epsilon$, use the fact that: $$\sqrt{a^2+b}\sim a + \frac{b}{2a} - \frac{b^2}{8a^3} + ...$$ So in our example, $$\sqrt{1732397825} \sim 41622.06575$$ Quite close to the true value of: $$\sqrt{1732397825}=41622.08338$$

-

You can use the Binomial Expansion to calculate square roots and other fractional powers. Since most square roots are irrational you get an infinite series, which you truncate to give you whatever level of accuracy/precision you require.

The Binomial Expansion only converges for |x| < 1 , so to get round this restriction take out a factor p where p is the largest perfect square ( or cube and so on... ) below s the number whose fractional index you are trying to evaluate.

For example for s = 5 the largest perfect square below it is 4 so

√5 = √4*√(1+1/4)      here x is 1/4 satisfying the condition for convergence.


The first few terms for √5 are:

√5 = √4[1 + 1/2*1/4 - 1/2*(1/4)²/2 + …… ]

√5 = 2(1   + 1/8     - 1/64        + …… )

-