# J.6 除法与余数 本节讨论浮点除法与余数。 ## 迭代除法 前面介绍过整数除法。将它改成浮点算法，与把整数乘法改成浮点乘法类似。由 $$ \frac{s_1\times2^{e_1}}{s_2\times2^{e_2}}=\frac{s_1}{s_2}\times2^{e_1-e_2} $$ 可知，只需用除法器计算 \(s_1/s_2\)，再乘以 \(2^{e_1-e_2}\)。图 J.2(b) 中的操作数装入位置与整数除法略有不同：\(s_2\) 装入 B，\(s_1\) 装入 P，A 无需保存操作数。随后使用整数除法算法，仅省去最初的一次左移，结果形式为 \(q_0.q_1\cdots\)。为舍入，再多算两个商位，分别作为保护位和舍入位，并用余数是否非零作为 sticky 位。保护位是必要的，因为第一个商位可能为 0；但分子分母均已规格化，商的最高两位不可能都为 0。这个算法每步产生一个商位。 另一种除法方法以二次而非线性速率收敛到商，J.10 节会介绍采用它的实际机器。下面先介绍两种主要迭代算法，再与直接除法算法比较。构造迭代算法的一般方法是 **Newton 迭代**，见图 J.13。先把问题变成求函数零点，再从某个零点估计出发，以该点的切线近似函数，用切线零点作为新估计。若 \(x_i\) 是估计值，切线方程为 $$ y-f(x_i)=f'(x_i)(x-x_i), $$ 其零点为 $$ x=x_{i+1}=x_i-\frac{f(x_i)}{f'(x_i)}.\tag{J.6.1} $$ ![用切线与横轴交点改进零点估计的 Newton 迭代](/figures/app-j/fig-j-13.png) {#fig-J.13 .wide .expand num="J.13" caption="用 Newton 迭代求零点。若 xi 是 f 的零点估计，xi+1 就是改进后的估计。求 f 在 xi 处的切线与 x 轴的交点，即得到 xi+1。"} 取 \(f(x)=x^{-1}-b\)，其零点为 \(1/b\)，所以 Newton 迭代给出从 \(b\) 求倒数的方法。因 \(f'(x)=-1/x^2\)，式 J.6.1 化为 $$ x_{i+1}=x_i-\frac{1/x_i-b}{-1/x_i^2} =x_i+x_i-x_i^2b=x_i(2-x_ib).\tag{J.6.2} $$ 因此，可如下计算 \(a/b\)： 1. 缩放 \(b\)，使 \(1\le b<2\)，查表取得 \(1/b\) 的近似值 \(x_0\)。 2. 迭代 \(x_{i+1}=x_i(2-x_ib)\)，直到 \(x_n\) 足够准确。 3. 计算 \(ax_n\)，并撤销第 1 步的缩放。 第 2 步需要多少次？精确到 \(p\) 位意味着相对误差的绝对值为 \(2^{-p}\)。简单代数计算表明，下一步的相对误差绝对值为 \(2^{-2p}\)，所以正确位数每步翻倍。Newton 迭代具有**自校正**性质：\(x_i\) 出现误差并不致命，它只是被当作 \(1/b\) 的一个估计，下一步会改进它，正确位数大致翻倍。舍入会造成误差；更重要的是，初始几次迭代本来就不期待很多正确位，因此可降低乘法精度，加快运算而不损害最终准确性。习题 J.20 讨论另一种应用。 > 译注：原页把下一步相对误差写为正的 \(2^{-2p}\)，但没有绝对值。令 \(\epsilon_i=bx_i-1\)，严格关系是 \(\epsilon_{i+1}=-\epsilon_i^2\)；正确位数翻倍的结论针对误差绝对值。 第二种方法常称为 **Goldschmidt 算法**。计算 \(a/b\) 时，让分子分母同时乘以使 \(rb\approx1\) 的数 \(r\)。具体令 \(x_0=a,y_0=b\)，每步计算 \(x_{i+1}=r_ix_i\)、\(y_{i+1}=r_iy_i\)，则 \(x_{i+1}/y_{i+1}=x_i/y_i=a/b\) 始终不变。选择 \(r_i\) 使 \(y_i\to1\)，便有 \(x_i\to a/b\)。 同样思想也能计算其他函数。例如，求 \(\sqrt a\)，令 \(x_0=y_0=a\)，每步计算 \(x_{i+1}=r_i^2x_i\)、\(y_{i+1}=r_iy_i\)，则 \(x_{i+1}/y_{i+1}^2=x_i/y_i^2=1/a\)。若选择 \(r_i\) 使 \(x_i\to1\)，就有 \(y_i\to\sqrt a\)。TI 8847 用这种技术求平方根。 回到除法。令 \(x_0=a,y_0=b\)，写 \(b=1-\delta\)，其中 \(|\delta|<1\)。若取 \(r_0=1+\delta\)，则 \(y_1=1-\delta^2\)；再取 \(r_1=1+\delta^2\)，则 \(y_2=1-\delta^4\)，依此类推。因为 \(|\delta|<1\)，所以 \(y_i\to1\)。相应有 $$ x_{i+1}=r_ix_i=(1+\delta^{2^i})x_i =\left[1+(1-b)^{2^i}\right]x_i, $$ 即 $$ x_{i+1}=a[1+(1-b)][1+(1-b)^2][1+(1-b)^4]\cdots[1+(1-b)^{2^i}].\tag{J.6.3} $$ 这个算法似乎有两个问题：\(b\) 离 1 较远，即 \(\delta\) 离 0 较远时，收敛慢；另外它不自校正，因为商是若干独立因子的乘积，其中一个因子的误差不会自行纠正。解决前者，可先查得 \(b\) 的近似倒数 \(b'\)，再对 \(ab'/bb'\) 运行算法；因为 \(bb'\approx1\)，收敛很快。 为应对后者，运算需多保留几位，以补偿舍入误差。不过，Goldschmidt 算法仍有较弱的自校正性质：\(r_i\) 的精确值并不重要。前几次迭代无需 \(1-\delta^{2^i}\) 的完整精度，可以把 \(r_i\) 取为 \(1+\delta^{2^i}\) 的截断值，以加快迭代而不影响收敛速度。此时 \(y_i\) 不再恰为 \(1-\delta^{2^i}\)，式 J.6.3 不能直接使用，但很容易组织成不依赖 \(r_i\) 精确值的算法： 1. 缩放 \(a,b\)，使 \(1\le b<2\)。 2. 查表取得 \(1/b\) 的近似值 \(b'\)。 3. 令 \(x_0=ab'\)、\(y_0=bb'\)。 4. 重复下列循环，直到 \(x_i\) 足够接近 \(a/b\)： $$ \begin{aligned} r&\approx2-y &&[y_i=1+\delta_i\ \Rightarrow\ r\approx1-\delta_i],\\ y&:=y\times r &&[y_{i+1}=y_i r\approx1-\delta_i^2],\\ x_{i+1}&:=x_i\times r &&[x_{i+1}=x_i r]. \end{aligned} $$ 两种迭代方法彼此相关。将 Newton 迭代展开，直接用 \(b\) 表示各项，而不递归引用 \(x_i\)，可得（见习题 J.22）： $$ x_{i+1}=x_0(2-x_0b)[1+(x_0b-1)^2][1+(x_0b-1)^4]\cdots[1+(x_0b-1)^{2^i}]. $$ 它与式 J.6.3 很相似。事实上，将式 J.6.3 中的 \(a,b\) 替换成 \(ax_0,bx_0\)，并取 \(a=1\)，两式完全相同。因此，若以无限精度迭代，两种方法会产生完全相同的 \(x_i\) 序列。 迭代的优点是无需专用除法硬件，可借用乘法器，只需增加控制逻辑；而且每步使正确位数翻倍，普通除法则每步只增加固定数量的位。 迭代求倒数有两个缺点。首先，IEEE 要求除法正确舍入，迭代却只能保证接近正确舍入值。Newton 方法先算 \(1/b\) 再乘 \(a\)，还存在额外问题：倒数正确舍入，并不保证商正确舍入。例如，十进制 \(p=2\)，\(a=13,b=51\)，则 \(a/b=.2549\ldots\) 应舍入为 .25；但 \(1/b=.0196\ldots\) 舍入为 .020，再乘 13 得 .26，末位相差 1。其次，迭代不给出余数；若借用浮点除法硬件做整数除法，这尤其麻烦，因为几乎所有高级语言都有余数运算。 传统经验认为，若想通过迭代得到正确舍入结果，应先将 \(1/b\) 算到略多于 \(2p\) 位，再将 \(a/b\) 算到略多于 \(2p\) 位，最后舍入到 \(p\) 位。但还有更快的方法，似乎最早实现于 TI 8847：仅增加约 6 位精度，算出初步商 \(q\)；再以同样额外 6 位精度比较 \(qb\) 与 \(a\)，就能迅速决定 \(q\) 是否正确舍入，或者最低位应加 1、减 1。习题 J.21 进一步讨论。 选择算法时，还要考虑除法与乘法的相对速度。除法更复杂，通常更慢。常见经验目标是使除法速度达到乘法的约三分之一。支持这一目标的理由之一，是一些实际程序，例如部分版本的 spice，除法与乘法次数之比为 1:3。另一个“3”来自平方根迭代：每次一次除法的算法可以换成每次三次乘法，见习题 J.20。 ## 浮点余数 对于非负整数，整数除法与余数满足 $$ a=(a\operatorname{DIV}b)b+a\operatorname{REM}b, \qquad0\le a\operatorname{REM}b 译注：这里余数界采用正除数 \(y\)；对一般非零实数除数，应写为 \(|y|/2\)。 若采用直接除法，求 REM 有自然办法：按计算 \(x/y\) 的过程进行。设 \(x=s_1 2^{e_1},y=s_2 2^{e_2}\)，在图 J.2(b) 中将 \(s_1\) 装入 P，\(s_2\) 装入 B。执行 \(e_1-e_2\) 次除法步骤后，P 保存 \(r=x-yn\)，满足 \(0\le r 译注：原页确实使用向下取整符号。但“每次至多 n 步、反复执行直到完成”通常还需计入不足 n 步的最后一次，不能把该式直接当作所有输入的精确调用次数。 REM 可用于三角函数。为简单起见，设十进制五位有效数字，计算 \(\sin x\)。若 \(x=7\)，可用 \(\pi\approx3.1416\) 缩减参数，算 \(\sin(7-2\times3.1416)=\sin(.7168)\)。但若 \(x=2.0\times10^5\)，则 \(2\times10^5/3.1416=63661.8\)，五位系统得到 63662；再乘 3.1416，精确积为 200000.5392，却舍入为 \(2.0000\times10^5\)，参数缩减竟得到 0，远非正确结果。问题是精度不足。 若有 REM，可算 \(2\times10^5\operatorname{REM}3.1416=-.53920\)，但仍不准确，因为 3.1416 只是 \(\pi\) 的近似值；实际 \(2\times10^5\operatorname{REM}\pi=-.071513\ldots\)。 传统上，计算大参数周期函数有两种做法：参数过大时报错，或者保存很多位 \(\pi\)，作精确参数缩减。REM 对这两者帮助都不大。某些数学库，例如 Berkeley UNIX 4.3bsd，使用第三种办法：把 \(\pi\) 舍入为最近浮点数，称为“机器 \(\pi\)”，记作 \(\pi'\)，然后用 \(x\operatorname{REM}\pi'\) 缩减参数。大 \(x\) 时它与 \(x\operatorname{REM}\pi\) 差别很大，因此这样计算并不能得到真实的 \(\sin x\)。不过，以这种方式计算各三角函数，可使 \(\sin^2x+\cos^2x=1\) 等熟悉恒等式在几个舍入误差范围内仍成立。因此，REM 配合机器 \(\pi\) 提供了一种简单办法：小参数时准确，大参数时仍可能有用。 参数缩减时，若 REM 还能返回 \(n\) 的低位，其中 \(x\operatorname{REM}y=x-ny\)，就很方便。实际实现往往用小于 \(2\pi\) 的量，例如 \(\pi/2\)，利用 \(\sin(x-\pi/2)=-\cos x\)、\(\sin(x-\pi)=-\sin x\) 等恒等式；这时需要 \(n\) 的低位来选择正确公式。