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 = D X_i - 1 \,,$$ then
 * $$X_i = {1 \over D} (1 + \epsilon_i) \,$$
 * $$\epsilon_{i+1} = - {\epsilon_i}^2 \,.$$