本文是对What Every Computer Scientist Should Know About Floating-Point Arithmetic的翻译。部分内容使用了ChatGPT先翻译然后校正而产生。原文的公式大部分是图片，本文都用mathjax重写了，方便阅读。作者加的注释都用【】括起来。
目录
摘要
浮点运算在许多人看来是一个神秘的课题。这相当令人惊讶，因为在计算机系统中，浮点运算无处不在。几乎每种编程语言都有浮点数数据类型；从个人电脑到超级计算机，都配备了浮点加速器；大多数编译器都会被要求编译浮点算法；而几乎每个操作系统都必须处理诸如溢出之类的浮点异常。本文介绍了浮点运算的一些方面，对计算机系统设计者产生直接影响。文章首先介绍了浮点表示(loating-point representation)和舍入误差(rounding error)的背景知识，然后讨论了IEEE浮点标准，并以大量示例结束，展示了计算机构建者如何更好地支持浮点运算。
简介
计算机系统的构建者经常需要关于浮点运算的信息。然而，关于这方面的详细信息确实非常有限。关于这个主题的少数几本书之一，帕特·斯特本茨（Pat Sterbenz）的《浮点运算(Floating-Point Computation)》早已绝版。本文是关于与系统构建直接相关的浮点运算（以下简称浮点）的教程。它由三个松散相关的部分组成。第一部分是“舍入误差”，讨论了在加法、减法、乘法和除法的基本操作中使用不同舍入策略的含义。它还包括关于两种舍入误差测量方法，即末尾单位(ulp)精度和相对误差(relative error)的背景信息。第二部分讨论了IEEE浮点标准，该标准正在迅速被商业硬件制造商接受。IEEE标准包括基本操作的舍入方法。对该标准的讨论基于“舍入误差”一节的材料。第三部分讨论了浮点与计算机系统各方面设计之间的联系。主题包括指令集设计、优化编译器和异常处理。
我尽量避免在谈论浮点时没有附上陈述为何正确的原因，特别是因为这些理由不涉及比基本微积分更复杂的内容。那些不是主要论点的解释已经被归纳到一个称为“细节”（The Details）的部分中，这样如果需要的话可以略过。特别是，许多定理的证明出现在这一部分。每个证明的结尾都标有符号z。当一个证明没有包括时，z会紧接着定理陈述之后出现。
舍入误差
将无限多的实数压缩到有限位数的比特中需要一种近似表示。虽然有无限多个整数，但在大多数程序中，整数计算的结果可以存储在32位中。相反，对于任何固定数量的比特，使用实数进行大多数计算会产生不能用这么多比特精确表示的量。因此，浮点计算的结果通常需要四舍五入以适应其有限的表示形式。这种舍入误差是浮点运算的特征性特点。”相对误差和末尾精度单位”一节描述了如何衡量这种舍入误差。
由于大多数浮点计算本身都存在舍入误差，如果基本算术运算引入比必要更多的舍入误差，这是否重要呢？这个问题贯穿于本节的主要主题。”保护位”一节讨论了保护位，这是一种在两个接近的数相减时减少误差的方法。IBM认为保护位非常重要，以至于在1968年将保护位添加到System/360架构的双精度格式中（单精度已经有保护位了），并为所有现有的计算机进行了改装。我们提供了两个示例来说明保护位的实用性。
IEEE标准不仅要求使用保护位，还提供了加法、减法、乘法、除法和平方根的算法，并要求实现产生与该算法相同的结果。因此，当程序从一台机器移植到另一台机器时，如果两台机器都支持IEEE标准，基本运算的结果将在每一位上都相同。这极大地简化了程序的移植过程。在”精确舍入操作”一节中还介绍了这一精确规范的其他用途。
浮点数格式
曾经有几种不同的实数表示方法被提出，但迄今为止最广泛使用的是浮点表示法。浮点表示法包括一个底数（通常假定为偶数）和一个精度 p。如果底数（基数）为 10，精度为 3，那么数字 0.1 的表示为 $1.00 × 10^{-1}$。如果底数为 2，精度为 24，那么十进制数 0.1 无法被精确表示，但可以近似表示为 $1.10011001100110011001101 × 2^{-4}$。
更一般的，一个浮点数可以表示为$\pm d.dd… d × \beta^e$，其中，d.dd…d 被称为尾数（significand），它有 p 位数字。更准确的说，$\pm d.dd… d × \beta^e$，表示：
\[\pm(d_0 + d_1\beta^{-1} + ... + d_{p-1}\beta^{-(p-1)})\beta^e, 0 \le d_i < \beta \;\;\;\;(1)\]术语”浮点数”指的是可以在讨论中精确表示的实数。与浮点表示法相关的另外两个参数是最大指数$e_{max}$和最小指数$e_{min}$。由于有p个可能的尾数（significand），以及$e_{max}-e_{min}+1$个可能的指数，因此可编码表示一个浮点数需要：$[log_2(e_{max}-e_{min}+1)]+[log_2\beta^p]+1$位，其中最后的 +1 是为了表示符号位。当前暂不需要关注精确的编码细节。
有两个原因导致一个实数可能无法被浮点数精确表示。最常见的情况可以用十进制数0.1来说明。尽管它有一个有限的十进制表示，在二进制中它有一个无限循环的表示。【每个有理数不管在什么进制里都是有限循环的，有限可以认为后面全是零的循环。比如0.1可以是0.100000…，也可以是0.0999999….，感兴趣的读者可以参考《什么是数学》。】因此，当$\beta$为2时，数字0.1严格位于两个浮点数之间，且不能被其中任何一个精确表示。较少见的情况是实数超出了范围，即其绝对值大于$\beta × \beta^e_{max}$或小于$1.0 × \beta^e_{min}$。本文大部分讨论了第一个原因导致的问题。然而，超出范围的数字将在“无穷大和非规格化数(Denormalized Number)”一节中进行讨论。
浮点表示法不一定是唯一的。例如，$0.01 × 10^1$ 和 $1.00 × 10^{-1}$ 都表示 0.1。如果最高位数字不为零（即方程（1）中的 $d_0 \ne 0$），则表示法被称为规格化。浮点数 $1.00 × 10^{-1}$是规格化的，而$0.01 × 10^1$不是。当$\beta=2$，精度$p=3$，$e_{min} = -1$，$e_{max} = 2$ 时，存在16个规格化浮点数，如下图所示:
粗体的标记(【上图可能看不清，其实就是数字1,2,4对于的那个竖线。】对应于尾数为 1.00 的数字。要求浮点表示法是规格化的使得表示法是唯一的。不幸的是，这个限制导致无法表示零！表示 0 的一种自然方式是 $1.0 × \beta^{e_{min}-1}$，因为这样可以保持非负实数的数值排序与其浮点表示的字典顺序对应的事实。当指数存储在 k 位字段中时，这意味着只有 $2^k - 1$ 个值可用作指数，因为必须保留一个值来表示 0。
请注意，浮点数中的乘号（×）是表示法的一部分，与浮点乘法运算不同。乘号（×）的含义应该根据上下文清楚。例如，表达式$(2.5 × 10^{-3})×(4.0 × 10^2)只涉及单一的浮点乘法运算。
相对误差和末尾单位
由于舍入误差是浮点计算中固有的，因此有必要有一种方法来衡量这种误差。考虑底数$\beta$为 10，精度p为3的浮点表示格式，该格式将贯穿本节。如果浮点计算的结果为 $3.12 × 10^{-2}$，而计算到无限精度时的答案是 0.0314，显然这个误差是末位(0.01)的2个单位(倍)。类似地，如果将实数 0.0314159 表示为 $3.14 × 10^{-2}$，那么它的误差为末位(0.01)的 0.159 个单位(倍)。通常，如果用浮点数 $d.d…d × \beta^e$ 来表示 z，那么其误差为 $|d.d…d - (z/\beta^e) |\beta^{p-1}$ 个末位单位。术语”ulps”将用作”末尾单位”的简称。如果计算的结果是离正确结果最近的浮点数，它的误差可能高达 0.5 ulp。衡量浮点数与其近似的实数之间差异的另一种方法是相对误差，它就是两个数字的差除以实数得到的结果。例如，用 3.14 × 100 来近似表示 3.14159 时产生的相对误差是 0.00159/3.14159 ≈ 0.0005。
【 补充一个例子说明一下。比如用$3.14 × 10^{-2}$来表示0.0314159，那么误差为$|3.14 - (0.0314159/10^{-2}) |10^{3-1}=0.159$。 或者用”人话来说，把误差的小数点先右移动e位再左移p-1位。比如上面的误差是0.0000159，右移-2位等价于左移2位就是0.00159，再左移2位就是0.159。 如果是0.0315呢，那么误差是0.0001，左移-2位然后左移2位就是1，所以用$3.14 × 10^{-2}$表示0.0315的误差是1个ulp。 】
要计算对应于0.5个末尾单位(ulp)的相对误差，可以观察到，当一个实数被最接近的可能浮点数 $d.dd…dd × \beta^e$ 近似时，误差可能高达 $0.00…0\beta’ × \beta^e$，其中$\beta’$是数字 $\beta/2$。【比如跟某个数最接近的是$3.14159 × 10^{-2}$，当这个数是$3.14159999…. × 10^{-2}$时误差最大，最大误差是$0.0000049999 × 10^{-2}$….=$0.000005 × 10^{-2}。$】浮点数的尾数中有 p 个单位，误差的尾数中有 p 个 0。这个误差为 $((\beta/2)\beta^{-p}) × \beta^e$。由于形如$d.dd…dd × \beta^e$的数字都具有相同的绝对误差，但其值范围在 $\beta^e$ 和$\beta × \beta^e$ 之间，相对误差范围在 $((\beta/2)\beta^{-p}) × \beta^e/\beta^e$和$((\beta/2)\beta^{-p}) × \beta^e/\beta^{e+1}$之间。也就是说，
\[\frac{1}{2} \beta^{-p} \le \frac{1}{2} ulp \le \frac{\beta}{2} \beta^{-p} \;\;\;\;(2)\]【 最后一句话值得解释一下，我当时想了很久才明白。所谓形如$d.dd…dd × \beta^e$指的是指数相同的浮点数，比如$1.00 × 10^{-2},1.01 × 10^{-2},…,9.99 × 10^{-2}$，这些数的最大误差都是$0.005 × 10^{-2}$，但是这些数的范围是$1.00 × 10^{-2} - 9.99 × 10^{-2}$。推广一下就是$\beta^e$ 和$\beta × \beta^e$之间，所以用最大绝对误差除以最大和最小的数就得到公式(2)。 另外我认为最大的数是不能取$\beta × \beta^e$的，比如上面的例子最大的是$9.99 × 10^{-2}$而不是$10 × 10^{-2}$，所以上面的第一个$\le $应该是严格小于。 】
特别地，对应于 0.5 ulp 的相对误差可以变化一个因子 。这个因子称为“wobble”。令$\epsilon = (\beta/2)\beta^{-p}$为（2）中上界的最大值，我们可以说，当一个实数四舍五入到最接近的浮点数时，相对误差始终受到$\epsilon$的限制，这被称为机器精度（machine epsilon）。
在上面的示例中，相对误差为 0.00159/3.14159 ≈ 0.0005。为了避免这样小的数值，通常将相对误差写成一个乘以$\epsilon$，而在这种情况下，$ \epsilon = (\beta/2)\beta^{-p}=5 × 10^{-3}= 0.005$。因此，相对误差可以表示为$(.00159/3.14159)/0.005 \epsilon \approx 0.1\epsilon$。
为了说明 ulps和相对误差之间的区别，考虑实数 $x = 12.35$。它近似为 \bar{x}= 1.24 × 10^1$。误差为 0.5 ulps，相对误差为 0.8。接下来考虑计算 $8x$。精确值为 $8x = 98.8$，而计算值为 $8\bar{x} = 9.92 × 10^1$。现在误差为 4.0 ulps，但相对误差仍然为 0.8。【相对误差从原来的$\frac{|\bar{x}-x|}{\bar{x}}$变成$\frac{|8\bar{x}-8x|}{8\bar{x}}，显然没有变化。】以 ulps 为单位测量的误差是相对误差的 8 倍，尽管相对误差相同。通常情况下，当底数为$\beta$时，以固定的相对误差表示的 ulps 可能会波动达到一个因子 $\beta$。相反，正如上述方程（2）所示，固定为 0.5 ulps 的误差会导致相对误差波动达到 $\beta$。
衡量舍入误差最自然的方式是以 ulps 为单位。例如，四舍五入到最接近的浮点数对应于小于或等于 0.5 ulp 的误差。然而，在分析各种公式引起的舍入误差时，相对误差是更好的衡量标准。这在”定理9”一节的分析中有很好的例证。因为$\epsilon$可能过高$\beta$倍的估计舍入到最近的浮点数，所以$\beta$较小的机器也会较好。【这也是二进制计算机的一个优势，因为最小的$\beta$就是2了！】
当只关心舍入误差的数量级时，ulps 和 $\epsilon$可以互换使用，因为它们之间最多相差一个因子$\beta$。例如，当浮点数的误差为 n ulps 时，这意味着受到影响的数字位数为 $log_\beta n$。如果计算中的相对误差为 $n\epsilon$，那么受影响的位数$\approx log_\beta n$。
【 我认为受影响的位数为$log_\beta n+1$。比如数字为$3.14 × 10^{-2}$，假设误差为1 ulps，也就是$3.13 × 10^{-2}$，那么受影响的位数就是最后一位。如果误差是10 ulps，那么就是$3.24× 10^{-2}$，影响2位。
相对误差如果是$n\epsilon$，因为$\epsilon=(\beta/2)\beta^{-p}$，所以$n\epsilon=n (\beta/2)\beta^{-p}$，而我们前面知道浮点数的范围是$\beta^p \le z \le \beta^{p+1}$，所以绝对误差的范围是$n (\beta/2)\beta^{-p} \times \beta^p$和$n (\beta/2)\beta^{-p} \times \beta^{p+1}$之间，也就是$n \beta/2$和 $n \beta^2/2$之间，那么绝对误差影响的位数就是把它取log，所以其位数是$log_\beta n +1 - log_\beta 2$和$log_\beta n + 2 - log_\beta 2$，这个值大约是$log_\beta n+1$(尤其是$\beta$较大比如10时)。 】
保护位
计算两个浮点数之间的差异的一种方法是精确计算差值，然后将其四舍五入到最接近的浮点数。如果操作数的大小差异很大，这种方法会非常耗时。假设 p = 3，计算 $2.15 × 10^{12} - 1.25 × 10^{-5}$ 可以表示为：
\[\begin{split} &x = 2.15 × 10^{12} \\ &y = .0000000000000000125 × 10^{12} \\ &x - y = 2.1499999999999999875 × 10^{12} \end{split}\]四舍五入结果为 $2.15 × 10^{12}$。与使用所有这些数字不同，浮点硬件通常在固定数量的数字上进行运算。假设保留的数字个数为 p，在较小的操作数向右移动时，数字会被简单丢弃（而不是四舍五入）。然后，$2.15 × 10^{12} - 1.25 × 10^{-5}$ 变为：
\[\begin{split} &x = 2.15 × 10^{12} \\ &y = 0.00 × 10^{12} \\ &x - y = 2.15 × 10^{12} \end{split}\]答案与如果将差计算为精确值然后四舍五入的结果完全相同。再举个例子：10.1 - 9.93。这将变为：
\[\begin{split} &x = 1.01 × 10^1 \\ &y = 0.99 × 10^1 \\ &x - y = .02 × 10^1 \end{split}\]正确答案是0.17，所以计算得到的差异高达30 ulps，每一位数字都是错误的！误差可以有多严重呢？
定理1
使用具有参数$\beta$和 p 的浮点格式，并使用 p 位数字计算差值时，结果的相对误差可能高达$\beta - 1$。
当x是1.00…0，y是$0.\rho\rho…\rho$时，其中$\rho=\beta-1$，相对误差最大。y的小数点后有p位，都是$\rho$。精确的差值是$x-y=\beta^{-p}$。但是当y只保留p位时，最后一位被忽略了，所以这样计算得到的差值是$\beta^{-p+1}$。所以误差是$|\beta^{-p}-\beta^{-p+1}|=\beta^{-p}(\beta-1)$，所以相对误差是$\beta^{-p}(\beta-1)/\beta^{-p}=\beta-1$，证毕。
当$\beta$为2时，误差可以和结果一样大；而当$\beta$为10时，误差可以是结果的9倍。换一种方式来看，当$\beta$为2时，公式(3)告诉我们，受影响的位数为$\log_2 1/\epsilon=log_22^p=p$。【相对误差为1，因此公式(3)的$n=1/\epsilon$。】也就是说，所有的位都错了。 为了解决这个问题，我们可以在运算的时候增加一位(这就是所谓的保护位)。也就是说，较小的数保留p+1位【较大数的第p位肯定是零】，运算结束后再舍入到p位。增加一个保护位后，前面的例子变成：
\[\begin{split} &x = 1.010 × 10^1 \\ &y = 0.993 × 10^1 \\ &x - y = .017 × 10^1 \end{split}\]上面的结果是完全准确的。【因为$.017 × 10^1$的规范化形式是$1.70 × 10^{-2}$，完全可以用规定的精度表示。】
但是在一个保护位的情况下，相对误差也可能大于$\epsilon$，比如110 - 8.59：
\[\begin{split} &x = 1.10 × 10^2 \\ &y = .085 × 10^2 \\ &x - y = 1.015 × 10^2 \end{split}\]舍入后变成102，而正确答案是101.41，相对误差.006，这比$\epsilon=\beta/2 × \beta^{-2}=5 × 10^{-2}=0.005$要大。不过一般情况下，相对误差只是比$\epsilon$大一点点。
未完待续！
- 显示Disqus评论(需要科学上网)