I am writing an algorithm that evaluates the square root of a positive real number $y$. To do this I am using the Newton-Raphton method to approximate the roots to $f(x)=x^2-y$. The $n^{th}$ iteration gives $$x_n=\frac{x_{n-1}^2+y}{2x_{n-1}}$$ as an approximation to $\sqrt{y}$. I found that starting with an initial guess $x_0=1$ works pretty well generally, so an answer to the question below that assumes $x_0=1$ is fine.

My question: is there an exact expression for the minimum $N$ of iterations needed to attain a given precision $p$ in the approximate solution $x_N$? In other words I'm looking for the smallest integer $N$ such that $$\left|\frac{x_N-\sqrt y}{\sqrt y}\right|<p.$$

I've thought about this for a while and played around with the expression for the errors $\epsilon_n = x_n - \sqrt y$ which can be shown to satisfy $\epsilon_{n+1}=\epsilon_n^2\,/\,2x_n$, but I can't find an answer. I've looked around on Google but I couldn't find an answer either.

Any pointers to a solution online or help would be much appreciated. A follow-up would of course be: can $x_0$ be optimised (while being a simple enough expression in terms of $y$) in order to minimise $N$?

  • $\begingroup$ I don't have an answer to your overall question, but here's a little bit of intuition: $x_n$ is the average of $x_{n-1}$ and $y/x_{n-1}$. One of these numbers must be less than $\sqrt{y}$, and one must be greater, so the average is a reasonable update. And of course, the closer $x_0$ is to $\sqrt{y}$, the better. Something like $x_0 = \lfloor \sqrt{y}\rfloor$ is a good guess (take the biggest perfect square less than $y$, and start with its square root). $\endgroup$ – BaronVT Nov 9 '13 at 16:58
  • $\begingroup$ Try analyzing $x_n^2 - y$ rather than $x_n - \sqrt{y}$. I imagine you probably want to analyze the roundoff error too, rather than assuming addition, multiplication and division are exact, which makes the problem more complicated. $\endgroup$ – Hurkyl May 2 '14 at 8:21

There is such a formula: consider

$$\frac{x_n+\sqrt y}{x_n-\sqrt y}=\frac{\frac{x_{n-1}^2+y}{2x_{n-1}}+\sqrt y}{\frac{x_{n-1}^2+y}{2x_{n-1}}-\sqrt y}=\frac{(x_{n-1}+\sqrt y)^2}{(x_{n-1}-\sqrt y)^2}=\left(\frac{x_{n-1}+\sqrt y}{x_{n-1}-\sqrt y}\right)^2.$$

By recurrence,

$$\frac{x_n+\sqrt y}{x_n-\sqrt y}=\left(\frac{x_{0}+\sqrt y}{x_{0}-\sqrt y}\right)^{2^n}.$$

If you want to achieve $2^{-b}$ relative accuracy, $x_n=(1+2^{-b})\sqrt y$,

$$2^n=\frac{\log_2\frac{(1+2^{-b})\sqrt y+\sqrt y}{(1+2^{-b})\sqrt y-\sqrt y}}{\log_2\left|\frac{x_{0}+\sqrt y}{x_{0}-\sqrt y}\right|},$$

$$n=\log_2\left(\log_2\frac{2+2^{-b}}{2^{-b}}\right)-\log_2\left(\log_2\left|\frac{x_{0}+\sqrt y}{x_{0}-\sqrt y}\right|\right).$$

The first term relates to the desired accuracy. The second is a penalty you pay for providing an inaccurate initial estimate.

If the floating-point representation of $y$ is available, a very good starting approximation is obtained by setting the mantissa to $1$ and halving the exponent (with rounding). This results in an estimate which is at worse a factor $\sqrt 2$ away from the true square root.

$$n=\log_2\left(\log_2\left(2^{b+1}+1\right)-\log_2\left(\log_2\frac{\sqrt 2+1}{\sqrt 2-1}\right)\right) \approx\log_2(b+1)-1.35.$$ In the case of single precision (23 bits mantissa), 4 iterations are always enough. For double precision (52 bits), 5 iterations.

On the opposite, if $1$ is used as a start and $y$ is much larger, $\log_2\left|\frac{1+\sqrt y}{1-\sqrt y}\right|$ is close to $\frac{2}{\ln(2)\sqrt y}$ and the formula degenerates to $$n\approx\log_2(b+1)+\log_2(\sqrt y)-1.53.$$ Quadratic convergence is lost as the second term is linear in the exponent of $y$.


If your starting value of $y$ is correctly "scaled" between $1$ and $100$, then a good initial guess is

$X_0 = 1,1545 + 0,11545*y$ (linear approx) or if you prefer a second order approx. you can use:

$X_0 = 1,0 + y*( 0,18 - 0,0009*y)$

Using $X_0$, you need only 3 iterations (Newton / Heron) to get a Relative Error less than 0,0001 % !

Now, here is a fast method to get an excellent $X_0$ (with less than 1% error):

Store a pre-calculated constants table containing all the square roots with rounded 4 decimals in the range $[1 ; 100]$. So $table(i)=sqrt(i)$.

Use Simple Linear Interpolation to get a good $X_0$. Now, the Maximum Relative Error is still greater than 1% ! ( 1.5% near $y=1.4142$).

You can simply divide this Maximum Relative Error (MRE) by a factor of 2 by doing this trick:

Simply replace the 1st entry ($\sqrt(1)=1.0000$ by $\sqrt(1)=1.0075$ ) AND the second entry ( $\sqrt(2)=1.4142$ by $\sqrt(2)=1.4248$ ). That's all. You have a MRE of 0,75%! So, the WORST case is a $X_0$ with only 0.75% in the WHOLE RANGE $[1 ; 100]$.

Unscale your number, and with only ONE Newton / Heron iteration get a 0.0028% MRE in any case, and with two iteration you have a 0.00000004% MRE!!

That's Fast, Simple and very Accurate. Have fun.

  • $\begingroup$ Please add Mathjax in your answers in future. $\endgroup$ – unseen_rider Mar 5 '17 at 10:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.