13.2  Reciprocal Refinement

A refinement $ y_{k+1}$ of a given reciprocal approximation $ y_k$ may be derived in two steps:

$\displaystyle e_k$ $\displaystyle =$ $\displaystyle \mathit{RNE}(1 - by_k, p)$  
$\displaystyle y_{k+1}$ $\displaystyle =$ $\displaystyle \mathit{RNE}(y_k + e_ky_k, p)$  

Alternatively, as illustrated by the double precision algorithm of Section 13.3, an approximation may be computed from the preceding two approximations as follows. This results in lower accuracy but allows the two steps to be executed in parallel:
$\displaystyle e_{k+1}$ $\displaystyle =$ $\displaystyle \mathit{RNE}(1 - by_{k+1}, p)$  
$\displaystyle y_{k+2}$ $\displaystyle =$ $\displaystyle \mathit{RNE}(y_k + e_ky_{k+1}, p)$  

The following lemma may be applied to either of these computations.

  (recip-refinement-1, recip-refinement-2) Assume that $ \vert 1 - by_1\vert \leq \epsilon_1$ and $ \vert 1 - by_2\vert \leq \epsilon_2$. Let $ e_1 = \mathit{RNE}(1 - by_1, p)$, and $ y_3 = \mathit{RNE}(y_3', p)$, where $ y_3' = y_1+e_1y_2$ and $ p > 0$. Let

$\displaystyle \epsilon'_3 = \epsilon_1(\epsilon_2 + 2^{-p}(1 + \epsilon_2)).
$

and

$\displaystyle \epsilon_3 = \epsilon'_3 + 2^{-p}(1 + \epsilon'_3).
$

Then (a) $ \vert 1 - by'_3\vert \leq \epsilon'_3$ and (b) $ \vert 1 - by_3\vert \leq \epsilon_3$.

PROOF: Let $ u_1 = 1 - by_1$ and $ u_2 = 1 - by_2$. Then $ \vert u_1\vert < \epsilon_1$, $ \vert u_2\vert < \epsilon_2$, and

$\displaystyle e_1 = (1 - by_1)(1 + v) = e_1(1 + v),
$

where $ \vert v\vert \leq 2^{-p}$. Thus,

$\displaystyle \vert 1 - by'_3\vert = \vert 1 - b(e_1y_2 + y_1)\vert = \vert(1 -...
...- u_1(1 + v)(1 - u_2)\vert = \vert u_1(u_2 + v(u_2 - 1))\vert \leq \epsilon'_3
$

and

$\displaystyle \vert 1 - by_3\vert \leq \vert 1 - by'_3\vert + b\vert y_3 - y'_3...
...^{-p}\vert by'_3\vert
\leq \epsilon'_3 + 2^{-p}(1 + \epsilon'_3)
= \epsilon_3. $

The inequality (b) of Lemma 13.2.1 provides a significantly reduced error bound for a refined reciprocal approximation $ y_3$ as long as the bounds $ \epsilon_1$ and $ \epsilon_2$ for the earlier approximations $ y_1$ and $ y_2$ are large in comparison to $ 2^{-p}$. To establish the bound $ 2^{-p}$ for the final approximation as required by Lemma 13.1.2, we shall use the inequality (a), pertaining to the corresponding unrounded value, in conjunction with the following additional lemma, which is a variation by Harrison [HAR00] of another result of Markstein [MAR90]. In practice, the application of this lemma involves explicitly checking a small number of excluded cases.

  (harrison-lemma) Let $ b$ be $ p$-exact. Assume that $ y = \mathit{RNE}(y', p)$ and $ \vert 1 - by'\vert \leq \epsilon' < 2^{-p-1}$ . Let $ d = \lceil 2^{2p}\epsilon'\rceil$. Then $ \vert 1 - by\vert < 2^{-p}$, with the possible exceptions $ b = 2 - 2^{1-p}k$, $ k = 1,\ldots,d$.

PROOF: If $ b = 1$, then $ \vert 1 - y'\vert < 2^{-p-1}$ implies $ y = 1$ and $ \vert 1 - by\vert = 0$. Thus we may assume $ b>1$, and therefore $ b \geq 1 + 2^{1-p}$. Consequently,

$\displaystyle y' \leq \frac{1}{b} (1 + 2^{-p-1}) < \frac{1 + 2^{-p-1}}{1 + 2^{1-p}} < 1,
$

and it follows that $ \vert y - y'\vert \leq 2^{expo(y')-p} \leq 2^{-p-1}$. Since

$\displaystyle \vert 1 - by'\vert \leq \epsilon' = 2^{-2p}2^{2p}\epsilon' \leq 2^{-2p}d,
$

and apart from the allowed exceptions, $ b < 2 - 2^{1-p}d$, we have

$\displaystyle \vert 1 - by\vert \leq \vert 1 - by'\vert + b \vert y - y'\vert
< 2^{-2p} d + (2 - 2^{1-p} d)2^{-p-1}
= 2^{-p}. $

David Russinoff 2017-08-01