附录D
What Every Computer Scientist
Should Know About Floating-Point
Arithmetic
注 – 本附录是对论文《What Every Computer Scientist Should Know About Floating-
Point Arithmetic》(作者:David Goldberg,发表于 1991 年 3 月号的《Computing
Surveys》)进行编辑之后的重印版本。版权所有 1991,Association for Computing
Machinery, Inc.,经许可重印。
D.1摘要
许多人认为浮点运算是一个深奥的主题。这相当令人吃惊,因为浮点在计算机系统中是普
遍存在的。几乎每种语言都有浮点数据类型;从 PC 到超级计算机都有浮点加速器;多
数编译器可随时进行编译浮点算法;而且实际上,每种操作系统都必须对浮点异常(如
溢出)作出响应。本文将为您提供一个教程,涉及的方面包含对计算机系统设计人员产生
直接影响的浮点运算信息。它首先介绍有关浮点表示和舍入误差的背景知识,然后讨论
IEEE 浮点标准,最后列举了许多示例来说明计算机生成器如何更好地支持浮点。
类别和主题描述符:(主要)C.0 [计算机系统组织]:概论—指令集设计;D.3.4 [程
序设计语言]:处理器—编译器,优化;G.1.0 [数值分析]:概论—计算机运算,错
误分析,数值算法(次要)
D.2.1 [软件工程]:要求/规范—语言;D.3.4 程序设计语言]:正式定义和理论—
语义;D.4.1 操作系统]:进程管理—同步。
一般术语:算法,设计,语言
其他关键字/词:非规格化数值,异常,浮点,浮点标准,渐进下溢,保护数位,NaN,
溢出,相对误差、舍入误差,舍入模式,ulp,下溢。
D-1
D.2简介
计算机系统的生成器经常需要有关浮点运算的信息。但是,有关这方面的详细信息来源非
常少。有关此主题的书目非常少,而且其中的一部《Floating-Point Computation》(作者:
Pat Sterbenz)现已绝版。本文提供的教程包含与系统构建直接相关的浮点运算(以下简
称为浮点)信息。它由三节组成(这三节并不完全相关)。第一节第 2 页的“舍入误差”
讨论对加、减、乘、除基本运算使用不同舍入策略的含义。它还包含有关衡量舍入误差的
两种方法 ulp 和相对误差的背景信息。第二节讨论 IEEE 浮点标准,该标准正被商业硬件
制造商迅速接受。IEEE 标准中包括基本运算的舍入方法。对标准的讨论借助了第 2 页的
“舍入误差”部分中的内容。第三节讨论浮点与计算机系统各个设计方面之间的关联。主
题包括指令集设计、优化编译器和异常处理。
作者已尽力避免在不给出正当理由的情况下来声明浮点,这主要是因为证明将产生较为
复杂的基本计算。对于那些不属于文章主旨的说明,已将其归纳到名为“详细信息”的
章节中,您可以视实际情况选择跳过此节。另外,此节还包含了许多定理的证明。每个证
明的结尾处均标记有符号❚。如果未提供证明,❚将紧跟在定理声明之后。
D.3舍入误差
将无穷多位的实数缩略表示为有限位数需要使用近似。即使存在无穷多位整数,多数程序
也可将整数计算的结果以 32 位进行存储。相反,如果指定一个任意的固定位数,那么多
数实数计算将无法以此指定位数精准表示实际数量。因此,通常必须对浮点计算的结果进
行舍入,以便与其有限表示相符。舍入误差是浮点计算所独有的特性。第 4 页的“相对
误差和 Ulp”章节说明如何衡量舍入误差。
既然多数浮点计算都具有舍入误差,那么如果基本算术运算产生的舍入误差比实际需要
大一些,这有没有关系?该问题是贯穿本章节的核心主题。第 5 页的“保护数位”章节
讨论保护数位,它是一种减少两个相近的数相减时所产生的误差的方法。IBM 认为保护
数位非常重要,因此在 1968 年它将保护数位添加到 System/360 架构的双精度格式(其
时单精度已具有保护数位),并更新了该领域中所有的现有计算机。下面两个示例说明了
保护数位的效用。
IEEE 标准不仅仅要求使用保护数位。它提供了用于加、减、乘、除和平方根的算法,并
要求实现产生与该算法相同的结果。因此,将程序从一台计算机移至另一台计算机时,如
浮点数的基数什么意思果这两台计算机均支持 IEEE 标准,那么基本运算的结果逐位相同。这大大简化了程序的
移植过程。在第 10 页的“精确舍入的运算”中介绍了此精确规范的其他用途。
D-2《数值计算指南》•2005 年 1 月
What Every Computer Scientist Should Know About Floating-Point Arithmetic D-3
D.3.1浮点格式
已经提议了几种不同的实数表示法,但是到目前为止使用最广的是浮点表示法。
1浮点表示法有一个基数 β(始终假定其为偶数)和一个精度 p 。如果 β = 10、p = 3,则将数 0.1表示为 1.00 × 10-1。如果 β = 2、p =24,则无法准确表示十进制数 0.1,但是它近似为
1.10011001100110011001101 × 2-4。
通常,将浮点数表示为 ± d.dd… d × βe ,其中 d.dd… d 称为有效数字2,它具有 p 个数字。更精确地说,± d 0 . d 1 d 2 … d p-1 × βe 表示以下数
。(1)术语浮点数用于表示一个实数,该实数可以未经全面讨论的格式准确表示。与浮点表示法相关联的其他两个参数是最大允许指数和最小允许指数,即 e max 和 e min 。由于存在 βp 个可能的有效数字,以及 e max – e min + 1 个可能的指数,因此浮点数可以按
位编码,其中最后的 +1 用于符号位。此时,精确编码并不重要。
有两个原因导致实数不能准确表示为浮点数。最常见的情况可以用十进制数 0.1 说明。虽然它具有有限的十进制表示,但是在二进制中它具有无限重复的表示。因此,当 β = 2时,数 0.1 介于两个浮点数之间,而这两个浮点数都不能准确地表示它。一种较不常见的情况是实数超出范围;也就是说,其绝对值大于 β × 或小于 1.0 × 。本文的大部分内容讨论第一种原因导致的问题。然而,超出范围的数将在第 20 页的“无穷”和第 22页的“反向规格化的数”章节中进行讨论。
浮点表示不一定是唯一的。例如,0.01×101 和 1.00 × 10-1 都表示 0.1。如果前导数字不
是零(在上面的等式 (1) 中,d 0 ≠ 0)
,那么该表示称为规格化。浮点数 1.00 × 10-1 是规格化的,而 0.01 × 101 则不是。当 β= 2、p = 3、e min = -1 且 e max = 2 时,有 16 个规格化浮点数,如图D-1 所示。粗体散列标记对应于其有效数字是 1.00 的数。要求浮点表示为规格化,则可以使该表示唯一。遗憾的是,此限制将无法表示零!表示 0 的一种自然方法是使用 1.0× ,因为 这样做保留了以下事实:非负实数的数值顺序对应于其浮
点表示的词汇顺序。
3将指数存储在 k 位字段中时,意味着只有 2k - 1 个值可用作指数,因为必须保留一个值来表示 0。
请注意,浮点数中的 × 是表示法的一部分,这与浮点乘法运算不同。× 符号的含义通过上下文来看应
该是明确的。例如,表达式 (2.5 × 10-3) × (4.0 × 102) 仅产生单个浮点乘法。
1.其他表示法的示例有浮动斜线和有符号对数 [Matula 和 Kornerup 1985;Swartzlander 和 Alexopoulos 1975]。
2.此术语由 Forsythe 和 Moler [1967] 提出,现已普遍替代了旧术语 mantissa 。
3.这假定采用通常的排列方式,即指数存储在有效数字的左侧。
d 0d 1β1–…d p 1–βp 1–()–+++⎝⎠⎛⎞β
e 0d i β<≤(),±log 2e max e min –1+()[]log 2βp ()[]1
++βe max βe min βe min 1–
D-4《数值计算指南》•2005 年 1 月
图D-1 β = 2、p = 3、 e min = -1、e max = 2 时规格化的数
D.3.2相对误差和 Ulp
因为舍入误差是浮点计算的固有特性,所以需要有一种方法来衡量此误差,这一点很重要。请考虑使用 β= 10 且 p = 3 的浮点格式(此格式在本节中广泛使用)。如果浮点计算的结果是 3.12 × 10-2,并且以无限精度计算的结果是 .0314,那么您可以清楚地看到最后一位存在 2 单位的误差。同样,如果将实数 .0314159 表示为 3.14 × 10-2,那么将在最后一位存在 .159 单位的误差。通常,如果浮点数 d.d …d × βe 用于表示 z ,那么将在最后一位存在 ⏐d.d …d − (z /βe )⏐βp-1 单位的误差。
4, 5 术语 ulp 将用作“最后一位上的单位数”的简写。如果某个计算的结果是最接近于正确结果的浮点数,那么它仍然可能存在 .5 ulp 的误差。衡量浮点数与它所近似的实数之间差值的另一种方法是相对误差,它只是用两数之差除以实数的商。例如,当 3.14159 近似为 3.14 × 100 时产生的相对误差是.00159/3.14159≈ .0005。
要计算对应于 .5 ulp 的相对误差,请注意当实数用可能最接近的浮点数 d.dd ...dd × βe 近似表示时,误差可达到 00β¢ × βe ,其中 β’ 是数字 β/2,在浮点数的有效数字中有 p
单位,在误差的有效数字中有 p 单位个 0。此误差是 ((β/2)β-p ) × βe 。
因为形式为 d.dd …dd × βe 的数都具有相同的绝对误差,但具有介于 βe 和 β × βe 之间的值,所以相对误差介于((β/2)β-p ) × βe /βe 和 ((β/2)β-p ) × βe /βe+1 之间。即,
(2)特别是,对应于 .5 ulp 的相对误差可能随因子 β 的不同而有所变化。此因子称为浮动系数。将 ε = (β/2)β-p 设置为上面 (2) 中的上限,表明:将一个实数舍入为最接近的浮点数
时,相对误差总是以 e (称为机器 ε)为界限。
在上例中,相对误差是 .00159/3.14159 ≈ .0005。为了避免此类较小数,通常将相对误差表示为一个因子乘以 ε 的形式,在此例中 ε = (β/2)β-p = 5(10)-3 = .005。因此,相对误差将表示为 (.00159/3.14159)/.005) ε ≈ 0.1ε。
以实数 x = 12.35 为例说明 ulp 与相对误差之间的差异。将其近似为 = 1.24 × 101。误差是 0.5 ulp ,相对误差是 0.8ε。接下来,将考虑有关数字 8的计算 。精确值是 8x =98.8,而计算值是 8 = 9.92 × 101。误差是 4.0 ulp ,相对误差仍然是 0.8ε。尽管相对误差保持不变,但以 ulp 衡量的误差却是原来的 9 倍。通常,当基数为 β 时,以 ulp 表示的固定相对误差可按最大因子 β 进行浮动。反之,如上述公式 (2) 所示,固定误差 .5 ulp 导致相对误差可按 β 进行浮动。
4.除非数 z 大于 +1 或小于 。否则在另行通知之前,将不会考虑以此方式超出范围的数。
5.假设 z ’ 是近似 z 的浮点数。那么,⏐d.d …d - (z /βe )⏐βp-1 将等价于 ⏐z ’-z ⏐/ulp(z ’)。衡量误差的更精确公式是 ⏐z ’-z ⏐/ulp(z )。– 编辑者
00.5
1234567βe max βe min 12--βp –12--ulp β2--βp
–≤≤x ˜x ˜x
˜
衡量舍入误差的最常用方法是以 ulp 表示。例如,舍入为最接近的浮点数相当于一个小于或等于 .5 ulp 的误差。但是,在分析由各种公式导致的舍入误差时,使用相对误差是一种较好的方法。在第 38 页的“证明”章节中的分析很好地说明了这一点。由于浮动因子β的关系,ε可能会过高估计舍入为最接近浮点数的效果,所以在使用较小β的计算机上,对公式的误差估计更为精确。
在仅关心舍入误差的大小顺序时,ulp 和ε可以互换使用,因为它们只是在因子β上有差
异。例如,在浮点数的误差为n ulp 时,这意味着受影响的位数是 log
β n。如果某个计算
的相对误差是nε,那么
受影响的数字是≈ log
β
n。(3) D.3.3保护数位
计算两个浮点数之差的一种方法是:精确计算二者之差,然后将其舍入为最接近的浮点
数。如果两个操作数的大小差别非常大,那么系统开销将也随之上升。假定p = 3,2.15
× 1012– 1.25 × 10-5将计算为
x = 2.15 × 1012
y = .0000000000000000125 × 1012
x-y = 2.1499999999999999875 × 1012
它舍入为 2.15 × 1012。浮点硬件通常按固定位数执行运算,而不是使用所有这些数字。假
定保留位数是p,并且在右移较小的操作数时,只是舍弃数字(与舍入相对)。那么,
2.15×1012–1.25×10-5将变为
x = 2.15 × 1012
y = 0.00 × 1012
x-y = 2.15 × 1012
结果是完全相同的,就好像先对差值进行精确计算,然后再进行舍入。请看另一个示例:
10.1 – 9.93。它将变成
x = 1.01 × 101
y = 0.99 × 101
x–y = .02 × 101
正确结果是 .17,因此计算差值的误差是 30 ulp,而且每一个数字均不正确:误差可以达
到多大程度?
D.3.3.1定理 1
使用带有参数β和 p的浮点格式,并使用p 数字计算差值,结果的相对误差可能与β– 1
一样大。
What Every Computer Scientist Should Know About Floating-Point Arithmetic D-5
版权声明:本站内容均来自互联网,仅供演示用,请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198,我们将在24小时内删除。
发表评论