Possible Double Rounding in Division Algorithms

Posted by Beetle B. on Thu 18 April 2019

This section deals with floating point \(a,b\) values, not necessarily between 1 and 2. Assume they are non-negative, though.

For this, let the significands be \(m_{a}\) and \(m_{b}\) respectively, and \(e_{a},e_{b}\) be the exponents. Use the self correcting iteration to calculate \(1/m_{a}\), and the other iteration to calculate \(m_{b}/m_{a}\). We’ll have \(q^{\ast}\) as a faithful approximation to \(m_{b}/m_{a}\), and \(y^{\ast}\) as a faithful approximation to \(1/m_{a}\) with relative error less than \(2^{-p}\).

Then apply the theorem to get \(r=\RN(m_{b}-m_{a}q^{\ast}\) and \(q'=\RN(q^{\ast}+ry^{\ast})\). Then \(q'=\RN(m_{b}/m_{a})\).

Note that \(q'\) is in the normal range. In fact, if \(a,b\) are normal numbers, it is in \((0.5,2)\) The smallest value even when you include denormals is \(1/(2^{p}-1)<2^{-p}\), which is still a normal number.

Let \(k=e_{b}-e_{a}\). If \(2^{k}q'\) is in the normal range, then \(\RN(b/a)=2^{k}q'\).

If it is larger than \(\Omega\), then return \(\infty\) for round to nearest (not sure if the author took into account the \(\frac{1}{2}\ulp(\Omega)\) aspect).

When \(2^{k}q'\) is a subnormal number, and is not a floating point number, this method may give the wrong value. The book talks about how to mitigate it - I didn’t go into the details.