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 迭代 — 图 J.13 用 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$ ：

缩放 $b$ ，使 $1\le b<2$ ，查表取得 $1/b$ 的近似值 $x_0$ 。
迭代 $x_{i+1}=x_i(2-x_ib)$ ，直到 $x_n$ 足够准确。
计算 $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$ 精确值的算法：

缩放 $a,b$ ，使 $1\le b<2$ 。
查表取得 $1/b$ 的近似值 $b'$ 。
令 $x_0=ab'$ 、 $y_0=bb'$ 。
重复下列循环，直到 $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<b.

类似地，可定义浮点余数为 $x=\operatorname{INT}(x/y)y+x\operatorname{REM}y$ 。怎样把 $x/y$ 转为整数？IEEE 余数采用就近取偶规则：选择满足 $|x/y-n|\le1/2$ 的整数 $n$ ；若有两个，取偶数。然后定义 $\operatorname{REM}=x-yn$ 。与整数余数范围不同，浮点余数满足 $|x\operatorname{REM}y|\le y/2$ 。这是精确定义，却不适合直接执行，因为 $n$ 可能极大；单精度中可达 $2^{127}/2^{-126}=2^{253}\approx10^{76}$ 。

译注：这里余数界采用正除数 $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<y$ 。IEEE 余数绝对值不超过 $y/2$ ，所以结果是 $r$ 或 $r-y$ 。只需保存最后一个商位，以处理恰好位于中点的情况。

但指数差可能很大，而浮点单元通常限制一条指令的最长执行时间，因此往往不能直接实现 REM。J.10 节的芯片都没有 REM，不过可像 Intel 8087 系列那样提供 remainder-step 指令：接收 $x,y$ ，不断执行除法步骤，直到 P 中得到余数，或已执行 $n$ 步； $n$ 是较小的常数，例如最高支持精度的除法步数。这样，软件可重复调用 REM-step，第一次用 $x$ 作分子，之后改用上一次的部分余数。原书给出的调用次数是 $\lfloor(e_1-e_2)/n\rfloor$ 。

译注：原页确实使用向下取整符号。但“每次至多 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$ 的低位来选择正确公式。