Number of iterations needed to attain a given precision $10^{-b}$ in Newton-Raphton method.

Question

I need to derive a formula to calculate the number of iterations needed to find $\sqrt{a}$ with a given accuracy of $10^{-8}$. The following is also given: $g(x)=\frac{1}{2}(x+\frac{a}{x})$, $a\in (1,100)$ and $x_0=10$. So

$$x_n=\frac{1}{2}\left(x_{n-1}+\frac{a}{x_{n-1}}\right)$$

How can I derive a formula to calculate the number of iterations needed? An answer to a similar question can be found here: Minimum number of iterations in Newton's method to find a square root. However the accuracy level there is $2^{-b}$ not $10^{-b}$. How should it be modified then? Also a big chunk of that proof isn't that clear to me. To be more precise (no pun intended) everything starting with

If you want to achieve $2^{−b}$ relative accura$\ldots$

is hazy (Why is $x_n=(1+2^{-b})\sqrt{y}$? Why is $2^n=\frac{\ldots}{\ldots}$?). Can anyone explain the linked answer or provide a better one? Any help would be much appreciated.

Answer 1

I thought the linked answer was quite elegant, but let's go back to basics. The relative error is $$\begin{align}e_n&=\left|\frac{\text{The value you've got}-\text{The exact value}}{\text{The exact value}}\right|\\ &=\left|\frac{x_n-\sqrt y}{\sqrt y}\right|=\frac{x_n-\sqrt y}{\sqrt y}\le10^{-b}\end{align}$$ Where I have dropped the absolute value bars because after the first iteration $x_n>\sqrt y$. This is equivalent to $x_n-\sqrt y\le10^{-b}\sqrt y$ or $x_n=(1+10^{-b})\sqrt y$.

Error refinement looks like this: given a relative error $e_n$ at step $n$ we see that $x_n=(1+e_n)\sqrt y$. Then $$\begin{align}x_{n+1}&=(1+e_{n+1})\sqrt y=\frac12\left(x_n+\frac y{x_n}\right)\\ &=\frac12\left((1+e_n)\sqrt y+\frac y{(1+e_n)\sqrt y}\right)\\ &=\frac12\left((1+e_n)\sqrt y+\sqrt y(1+e_n)^{-1}\right)\\ &\approx\frac{\sqrt y}2\left(1+e_n+1-e_n+e_n^2\right)=\sqrt y\left(1+\frac12e_n^2\right)\end{align}$$ It follows that $e_{n+1}\approx\frac12e_n^2$. Taking natural logarithms, $\ln e_{n+1}=2\ln e_n-\ln2$. If we let $Z_n=\ln e_n$ we have a linear difference equation for $Z_n$: $$Z_{n+1}=2Z_n-\ln2$$ We first solve the homogeneous equation: $$Z_{n+1,h}=2Z_{n,h}$$ By assuming a solution of the form $Z_{n,h}=Cr^n$. This leads to $Cr^{n+1}=2Cr^n$ or $r=2$, so $Z_{n,h}=C2^n$. Then we find a particular solution to the equation; we will guess $Z_{n,p}=A$, a constant. Then $$Z_{n+1,p}=A=2Z_{n,p}-\ln2=2A-\ln2$$ With solution $Z_{n,p}=A=\ln2$. So now we put together $$\ln e_n=Z_n=Z_{n,h}+Z_{n,p}=C2^n+\ln2$$ Applying the initial conditions $\ln e_0=C+\ln2$ we find that $C=\ln\frac{e_0}2$ so $$\ln e_n=\left(\ln\frac{e_0}2\right)2^n+\ln2$$ And finally $$e_n=2\left(\frac{e_0}2\right)^{2^n}\le10^{-8}$$ It's possible to get an initial estimate with relative error $e_0=0.03$ thus $2(0.015)^{2^n}\le10^{-8}$ so $(0.015)^{2^n}\le5\times10^{-9}$ or $2^n\ln0.015\le\ln5\times10^{-9}$ or $2^n\ge\frac{\ln5\times10^{-9}}{\ln0.015}$ or $$n\ge\frac{\ln\left(\frac{\ln5\times10^{-9}}{\ln0.015}\right)}{\ln2}\approx2.19$$ So we need $3$ iterations to get to relative error of $10^{-8}$. I don't have time right now to explain my original approximation, here is a test program:

module sqrtmod
   use ISO_FORTRAN_ENV
   implicit  none
   private
   public mysqrt
   contains
      function mysqrt(x) bind(C,name='mysqrt')
         real(REAL64) mysqrt
         real(REAL64),value :: x
         integer(INT64), parameter :: magic = int(Z'3feed9eba10e6360',INT64)
         integer i
!! Initial approximation I didn't have time to explain !!
         mysqrt = transfer(shiftr(magic+transfer(x,magic),1),mysqrt)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         do i = 1, 3
            mysqrt = 0.5_REAL64*(mysqrt+x/mysqrt)
         end do
      end function mysqrt
end module sqrtmod

program sqrt1
   use ISO_FORTRAN_ENV
   use sqrtmod
   implicit  none
   integer(INT64) magic
   real(REAL64) x, y, maxerr
   integer(INT64) i, imax
   imax = 10000000
    maxerr = 0
   do i = 0, imax
      x = 1+i*4.0_REAL64/imax
      y = mysqrt(x)
      maxerr = max(maxerr,(y-sqrt(x))/sqrt(x))
   end do
   write(*,*) 'maxerr = ',maxerr
end program sqrt1

EDIT: A little free time. To get the initial approximation we observe that $$\sqrt y=y^{1/2}=2^{\frac12\log_2y}$$ So if we could somehow take the base $2$ logarithm of $y$, divide it by $2$, and then take $2$ to that power, we would have a good approximation to, in fact exactly, $\sqrt y$. Well, as an $\text{IEEE-754}$ double precision floating point number, the base $2$ logarithm is available as bits $52:62$ of $y$. To divide by $2$ we just shift the bits of the representation right by $1$. Then raising $2$ to that power means just interpreting the result as a floating point number once again.

Now, it's a little more complicated than that because the exponent has a bias and there is an implicit leading bit in the mantissa and shifting the mantissa right by $1$ doesn't take its square root! So we have to add a magic number to the integer that represents $y$ as a floating point number to even get in the ballpark but the initial value $x_0$ that we end up with is within about $3\%$ of the exact value.

Let's show a little example of how this works. Start with $y=\frac{\pi^2}8$ $$\begin{align} y &= 1.2337005501361697\\ y &= \mathtt{\text{ 3FF3BD3CC9BE45DE}}\\ \text{magic} &= \mathtt{\text{ 3FEED9EBA10E6360}}\\ y+\text{magic} &= \mathtt{\text{ 7FE297286ACCA93E}}\\ (y+\text{magic})/2 &= \mathtt{\text{ 3FF14B943566549F}}\\ x_0 &= 1.0809518896033194\\ \sqrt y &= 1.1107207345395915\\ \text{% error} &= 2.68\\ \end{align}$$