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:
- 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$.
- 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:
- 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,
floororceilfunctions, or any other basic tools. - 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.
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) + 0x3F7A3BEAif I scaled the magic number correctly.