The general form of Newton's Method for zero-finding is
\begin{equation*}
x_{k+1} := x_k - \frac{f(x_k)}{f'(x_k)}
\end{equation*}
The Reciprocal Square Root Newton algorithm is obtained by setting $f(x) = \frac{1}{x^2-a}$ and we get
\begin{equation*}
x_{k+1} := x_k - \frac{f(x_k)}{f'(x_k)} = x_k - \frac{1/x_k^2-a}{-2/x_k^3}=x_k +
\frac{1}{2}x_k(1-ax^2_k).
\end{equation*}
The sequence $\{x_0,x_1, \ldots, x_k, \ldots\}$ converges quadratically, for all $x_0 >0$, to a fixed point $x$. That is
\begin{equation*}
x = x + \frac{1}{2}x(1-ax^2),
\end{equation*}
which has two solutions, $x = 0$, and $x = 1/\sqrt{a}$.
We compute $x_{k+1}$ in two steps:
- $e_k := 1-a\star x_k\star x_k $
- $x_{k+1} := x_k + e_k\star x_k/2$
Here, $e_k$ is the relative error in $x^2_k,$ i.e.,
$(1/a-x^2_k)/(1/a)=(1-ax^2_k),$ and $e_kx_k/2$ is the correction term.
The most important aspect of this algorithm is that no divisions are
required. The division by 2 can be done by bit-shifting. Calculating $x_{k+1}$ requires 3 Mults
and 2 Adds
. One final multiplication, $a\star x_k$, is done to calculate $\sqrt{a}$.
If a suitable starting value for $x_0$ is chosen, then this algorithm converges to full IEEE double precision in about 5 iterations. An algorithm similar to this is used by Intel in their Itanium processor, whose basic arithmetic operation is the fused multiply-add (FMA) instruction, $z = ax+y$. All other operations are performed using this FMA instruction. For historical reasons, Intel shies away from long division.
See Peter Markstein, IA-64 and Elementary Functions, Prentice-Hall, 2000
Additional Information
The Reciprocal Newton algorithm is obtained by setting $f(x) = \frac{1}{x-a}$ and we get
\begin{equation*}
x_{k+1} := x_k - \frac{f(x_k)}{f'(x_k)} = x_k - \frac{1/x_k-a}{-1/x_k^2}=x_k +
x_k(1-ax_k).
\end{equation*}
We compute $x_{k+1}$ in two steps:
- $e_k := 1-a\star x_k $
- $x_{k+1} := x_k + e_k\star x_k/2$
Here, $e_k$ is the relative error in $x^2_k,$ i.e.,
$(1/a-x_k)/(1/a)=(1-ax_k),$ and $e_kx_k/2$ is the correction term. Each of these two steps is calculated by a single FMA instruction. The division, $z = y/x$, is calculated as $y\star\frac{1}{x}$. Note how similar this algorithm is to the reciprocal square root algorithm.
Choosing a starting value $x_0$ for the Reciprocal Algorithm is a delicate matter. It can be shown that the domain of convergence is $ 0 < x_0 < \frac{2}{a}$, but this begs the question: what is $\frac{2}{a}$? This is getting into deep water, so I'll leave the answer to Peter Markstein above.