User:PeterDixon77/sandbox

In order to apply Newton's method to find the reciprocal of $D$, it is necessary to find a function $f(X)$ which has a zero at $X=1/D$. The obvious such function is $f(X)=DX-1$, but the Newton–Raphson iteration for this is unhelpful since it cannot be computed without already knowing the reciprocal of $D$. Moreover multiple iterations for refining reciprocal are not possible since higher order derivatives do not exist for $f(X)$. A function which does work is $f(X)=1/X-D$, for which the Newton–Raphson iteration gives

$X_{i+1}=X_{i}-{f(X_{i}) \over f'(X_{i})}=X_{i}-{1/X_{i}-D \over -1/X_{i}^{2}}=X_{i}+X_{i}(1-DX_{i})=X_{i}(2-DX_{i})$

which can be calculated from $X_{i}$ using only multiplication and subtraction, or using two fused multiply–adds.

From a computation point of view the expressions $X_{i+1}=X_{i}+X_{i}(1-DX_{i})$ and $X_{i+1}=X_{i}(2-DX_{i})$ are not equivalent. To obtain a result with a precision of n bits while making use of the second expression one must compute the product between $X_{i}$ and $(2-DX_{i})$ with double the required precision (2n bits). In contrast the product between $X_{i}$ and $(1-DX_{i})$ need only be computed with a precision of n bits.

If the error is defined as $\epsilon _{i}=DX_{i}-1\,,$ then

$X_{i}={1 \over D}(1+\epsilon _{i})\,$

$\epsilon _{i+1}=-{\epsilon _{i}}^{2}\,.$