Algorithm to estimate $\sqrt{a}$ for given $a\in\mathbb{Z}^+$

71 Views Asked by At

Given positive integer $a$, I want to compute its square root (without using sqrt function). I am using Newton Raphson method to do it. Our function is $f(x)=x^2-a$, where $x\in\mathbb{R}$.

$$x_{n+1}=x_{n}-\frac{f(x_{n})}{f'(x_{n})}.$$

I know that Newton Raphson converges fast and for this example the value of the initial guess is not that important. But anyway, I start wondering about finding the best initial guess $x_0$ using simple math tools (no radicals, no differentiation).

Questions:

  1. Given any number $a$, what is the best guess of $x_0$, such that $x_0$ is an integer and $x_0^2-a$ is as close as possible to $0$.
  2. Given any number $a$, what is the best guess of $x_0$, such that $x_0$ is a real and $x_0^2-a$ is as close as possible to $0$.

There are few things I want to mention:

  1. Remember that we don't know how to compute the value of $\sqrt{a}$, so algorithm should reflect it. We can use basic math tools as $+,-,\times$ and $/$, or tools like counting number of digits, floor or ceil functions, or any other basic tools.
  2. We cannot assume $a=3$ for example, and start estimating $x_0$. For any given $a$, I want to follow the same algorithm for the best guess.
1

There are 1 best solutions below

3
On BEST ANSWER

If you want to be pedantic, then any option you give for the "best" first guess can instantly be improved by taking that and iterating it through Newton-Raphson. So in other words, if you tell me you're starting with $m$, I'll start with $\frac{m^2 + a}{2m}$ and be one step ahead of you. Or I could nest a few iterations together and get something even better.

And since Newton-Rahpson converges so quickly (roughly doubling the number of correct digits per iteration), it's very hard to find a good starting value that isn't more expensive in terms of operations than just choosing a mediocre starting point and jumping right into the algorithm. So on that basis, something as simple as $x_0 = a$ is perfectly fine.

For numbers that are going to be stored in floating-point format, you could possibly use a similar hack to Quake's fast inverse square root algorithm, since that uses very low-level binary operations to get its first approximation. The equivalent expression for just taking the square root would be $I_y \approx I_x + L(B - \sigma)$, using the notation from the Wikipedia article, which in code form would be i = (i >> 1) + 0x3F7A3BEA if I scaled the magic number correctly.