Oracle® Solaris Studio 12.4:数值计算指南

退出打印视图

更新时间: 2015 年 1 月
 
 

What Every Computer Scientist Should Know About Floating-Point Arithmetic》补充资料

在阅读本数值计算指南时,每个用户都可以在 David Goldberg 于 1991 年 3 月发表在《Computing Surveys》上的论文《What Every Computer Scientist Should Know About Floating-Point Arithmetic》中找到有用的信息(请参见 http://dl.acm.org/citation.cfm?id=103163)。

论文摘要如下:

许多人将浮点运算看作非常深奥的话题。这相当让人吃惊,因为浮点在计算机系统中无所不在:几乎每种语言都有浮点数据类型;从 PC 到超级计算机在内的计算机都具有浮点加速器;绝大多数编译器都会被不断调用来编译浮点算法;几乎所有操作系统都必须响应诸如溢出等浮点异常。本文提供了教程,介绍浮点对计算机系统的设计人员具有直接影响的各个方面。文中首先介绍了浮点表示形式和舍入误差的背景信息,然后对 IEEE 浮点标准进行了讨论,并通过示例来总结计算机系统开发人员如何更好地支持浮点运算。

此补充资料不是已发布的 Goldberg 论文的一部分。增加此内容的目的是说明一些要点,并更正读者在阅读论文时对 IEEE 标准可能产生的一些不正确概念。本材料并非由 David Goldberg 撰写,但在此处使用已得到他的授权。标准 754-1985 和 854-1987 已由 754-2008 取代,后者指定了二进制和十进制浮点运算。这不影响 Goldberg 的论文。

此补充资料专门讨论了 IEEE 754 实现之间的差别。本主题包含以下子主题:

D.1 IEEE 754 实现之间的差别

Goldberg 论文指出必须谨慎实现浮点运算,因为程序员需要依靠其属性来确保程序的正确性和精确度。特别是,IEEE 标准要求谨慎实现,只有在符合标准的系统上才能够编写能够正常工作并提供准确结果的有用程序。读者可能会得出结论,这样的程序应该可以在所有 IEEE 系统上移植。实际上,如果“程序在两台均支持 IEEE 运算的计算机之间移动时,只要出现任何中间结果差异,一定是因为软件错误而不是算法的不同所造成的。”这种说法是正确的,那么编写可移植软件会容易得多。

不幸的是,IEEE 标准不确保相同的程序在所有符合标准的系统上提供完全一致的结果。在不同的系统上,由于各种原因,大部分程序实际上会生成不同的结果。例如,许多程序使用系统库提供的初等函数,而在标准中并未完全指定这些函数。1985 标准未完全指定数值在十进制和二进制格式之间的转换,也完全没有指定超越函数。

许多程序员可能未认识到,即使程序只使用 IEEE 标准中规定的数字格式和运算,在不同的系统上也会得到不同的结果。实际上,标准的作者的目的是允许不同的实现获得不同结果。在 IEEE 754 标准中定义术语目标时,其意图非常明显:“目标可以由用户显式指定,也可以由系统隐式指定(例如,子表达式中的中间结果或过程的参数)。一些语言将中间计算的结果放在用户无法控制的目标中。尽管如此,此标准按照该目标的格式和操作数的值来定义运算的结果。”(IEEE 754-1985,第 7 页)换而言之,IEEE 标准要求按照其中将放入结果的目标的精度,正确地舍入每个结果,但是该标准不要求目标的精度由用户程序确定。因此,不同系统可能会将其结果使用不同精度提供给目标,导致相同的程序生成不同结果(有时是动态生成),即使这些系统均符合标准也是这样。

引用论文中的多个示例运用了浮点运算舍入方法的一些知识。要想依靠这样的一些示例,程序员必须能够预测解释程序的方式,特别是在 IEEE 系统上,每个算术运算目标的精度是多少。IEEE 标准中关于目标定义的漏洞会影响到程序员理解如何解释程序的能力。因此,Goldberg 论文中提供了多个示例,这些示例在高级语言中作为明确可移植的程序来实现时,可能在 IEEE 系统上无法正确运行,因为这些系统通常使用与程序员预期所不同的精度向目标提供结果。另一个示例也许能运行,但也证明,程序的正常运行可能超出了程序员的一般能力范围。

在本补充资料中,IEEE 754 运算的当前实现基于通常所用的目标格式精度进行了分类。在回顾论文中的一些示例时,表明如果所提供结果的精度比程序预期的精度更高,会导致程序计算得出错误的结果,即使可以证明在使用预期精度时是正确的。论文中修订了这些证据之一,说明即使这并不表明程序是错误的,也应该掌握一些处理意外精度所需的知识技能。这些示例表明,尽管 IEEE 标准已有规定,其中允许的不同实现之间的差异会妨碍您编写可移植、高效且可准确预测结果的数值软件。要开发这样的软件,首先必须创建编程语言和环境,限制 IEEE 标准允许的可变性,从而使您能够表达出程序所需的浮点语义。1985 标准的 2008 版本针对编程语言做出了这样的推荐。

D.1.1 当前 IEEE 754 实现

当前 IEEE 754 运算的实现可以分为两组,按照它们在硬件中支持不同浮点格式的程度进行区分。扩展精度系统已使用 Intel x86 处理器系列和 –xarch=386 选项编译得到例证,为扩展双精度格式提供了完整支持,但对单精度和双精度只提供了部分支持:它们提供了指令,用于以单精度或双精度加载或存储数据,将其实时在扩展双精度格式之间来回转换,并提供特殊模式(而非缺省模式),其中数学运算的结果将舍入为单精度或双精度,即使将结果以扩展双精度格式保存在寄存器中也是如此。在这些模式中,Motorola 68000 系列处理器将结果舍入到单精度或双精度格式的精度和范围。在 Intel x86 和兼容处理器中,将结果舍入到单精度或双精度格式的精度,但保留与扩展双精度格式相同的范围。单/双精度系统,包括大多数 RISC 处理器,提供了对单精度和双精度格式的完整支持,但不支持符合 IEEE 标准的扩展双精度格式。x86 SSE2 扩展提供了 SSE2 寄存器中的单/双精度系统,仍支持扩展精度寄存器中的扩展精度。

要查看计算在扩展精度系统上相比单/双精度系统可能会有什么不同的行为,请考虑 Goldberg 论文中示例 "Systems Aspects" 的 C 版本,如下所示:

int main() {
    double  q;
 
    q = 3.0/7.0;
    if (q == 3.0/7.0) printf("Equal\n");
    else printf("Not Equal\n");
    return 0;
}

常量 3.0 和 7.0 解释为双精度浮点数,表达式 3.0/7.0 继承 double 数据类型。在单/双精度系统上,表达式将以双精度求值,因为这是可以使用的最有效格式。因此,将向 q 分配正确舍入到双精度的值 3.0/7.0。在下一行中,表达式 3.0/7.0 将以双精度格式再次求值,结果将等于刚赋值给 q 的值,因此程序将按预期输出 "Equal"(等于)。

在扩展精度系统上,即使表达式 3.0/7.0 具有类型 double,商将以扩展双精度格式在寄存器中计算,因此在缺省模式下,它将舍入到扩展双精度。但是,将得到的值赋值给变量 q 时,它可能会存储在内存中,由于 q 声明为 double,该值将舍入为双精度。在下一行中,表达式 3.0/7.0 可使用扩展精度再次求值,得到的结果不同于存储在 q 中的双精度值,这使得程序输出 "Not equal"(不等于)。

也有可能出现其他结果:编译器可以决定存储,从而舍入第二行中表达式 3.0/7.0 的值,然后再将其与 q 比较,或者可以将 q 以扩展精度格式保存在寄存器中而不存储。优化编译器可能会在编译时对表达式 3.0/7.0 求值,可能采用双精度,也可能采用扩展双精度。使用一个 x86 编译器,程序在采用优化编译的情况下输出 "Equal"(等于),在针对调试进行编译时输出 "Not Equal"(不等于)。最后,一些面向扩展精度系统的编译器会自动更改舍入精度模式,导致在寄存器中生成结果的运算将这些结果舍入到单精度或双精度,尽管可能具有更宽的范围。因此,在这些系统上,只是阅读程序的源代码并运用对 IEEE 754 运算的基本理解是无法预测其行为的。无法提供符合 IEEE 754 要求的环境的原因不能归结于硬件或编译器。正如对硬件要求的那样,硬件向每个目标提供正确舍入的结果,而编译器则将一些用户无法控制的中间结果分配给目标,这也是系统所允许的。

这一节的剩余部分将介绍使用 –xarch=386 编译的 x86 行为,这是 Oracle Solaris Studio 12 及早期发行版中的缺省行为。在 Oracle Solaris Studio 12.4 中,缺省值为 –xarch=sse2

D.1.2 在扩展精度系统上执行计算的陷阱

传统观点认为,相比单/双精度系统,扩展精度系统如果不能做到更精确,至少应该生成相同精确度的结果,因为扩展精度系统至少能够提供相同精度,通常更高。以下章节中讨论的普通示例以及基于这些示例的更精细的程序将表明,这种想法完全是一厢情愿:一些明确可移植的程序,确实可以跨单/双精度系统进行移植,但在扩展精度系统上会提供不正确的结果,恰恰是因为编译器和硬件有时提供的精度会超过程序预期的精度。

目前的编程语言使得程序难以指定预期的精度。如 Goldberg 论文中的 "Languages and Compilers" 一节所述,许多编程语言没有规定 10.0*x 这样的表达式在相同上下文中每次出现时都应该得到相同的值。在 IEEE 标准之前,一些语言(例如 Ada)由于不同运算之间的差别会受到这些方面的影响。近来,诸如 ANSI C 等语言已受到符合标准的扩展精度系统的影响。实际上,ANSI C 标准明确允许编译器对浮点表达式的求值结果使用高于通常与结果类型关联的精度。因此,表达式 10.0*x 的值可能会由于各种因素而变化:表达式是立即赋值给变量还是作为较大表达式中的子表达式;表达式是否参与比较;表达式是否作为参数传递到函数,如果是,参数是按值传递还是按引用传递;当前精度模式;编译程序的优化级别;编译器在编译程序时使用的精度模式和表达式求值方法等等。

表达式求值时的异常行为不能完全归咎于语言标准。只要表达式在扩展精度寄存器中求值,扩展精度系统的运行效率最高,即使需要存储的值以所需的最小精度存储。约束语言以便要求 10.0*x 在任何位置求值时得到相同值会对这些系统的性能造成不利影响。不幸的是,允许这些系统在语法等同的上下文中以不同方式对 10.0*x 求值时,自身会对准确数值软件的程序员造成不利影响,使他们无法依靠程序的语法来表达出所希望的语义。

以下示例探讨了真实程序是否能依赖于给定表达式始终可以计算得到相同值的假设。请回顾一下 Goldberg 论文中有关计算 ln(1 + x) 的定理 4,这是用 Fortran 编写的:

real function log1p(x)
real x
if (1.0 + x .eq. 1.0) then
   log1p = x
else
   log1p = log(1.0 + x) * x / ((1.0 + x) - 1.0)
endif
return

在扩展精度系统上,编译器可能会在扩展精度的第三行对表达式 1.0 + x 求值,并将结果与 1.0 进行比较。但是,将相同的表达式传递到第六行中的对数函数时,编译器可能会将其值存储在内存中,舍入到单精度。因此,如果 x 不是太小(在扩展精度下 1.0 + x 才会舍入到 1.0),而是足够小(在单精度下 1.0 + x 可以舍入到 1.0),则 log1p(x) 返回的值将为零而不是 x,相对误差为 1,而不是大于 5ε。与此类似,假定第六行表达式的剩余部分(包括重复出现的子表达式 1.0 + x 在内)在扩展精度下求值。在这种情况下,如果 x 很小,但又不足够小,1.0 + x 在单精度下舍入到 1.0,则 log1p(x) 返回的值会超过正确值达到 x,相对误差仍接近 1。如果需要明确的示例,可以让 x 为 2-24 + 2-47,因此 x 是最小的单精度数,这样 1.0 + x 向上舍入到下一个较大数 1 + 2-23。然后,log(1.0 + x) 约为 2-23。由于表达式中第六行的分母在扩展精度下求值,将对其精确求值并提到 x,因此 log1p(x) 返回的值大约为 2-23,这接近精确值的两倍。

至少一个编译器上会实际出现这种情况。使用面向 x86 系统的 Sun WorkShop Compilers 4.2.1 Fortran 77 编译器和 -O 优化标志编译上述代码时,生成的代码精确按照所述方式计算 1.0 + x。因此,函数将 log1p(1.0e-10) 计算为 0,将 log1p(5.97e-8) 计算为 1.19209E-07

要使定理 4 的算法正确运行,表达式 1.0 + x 必须在每次出现时以相同方法求值。仅当 1.0 + x 在一个实例中求值为扩展双精度,在另一个实例中求值为单精度或双精度时,该算法在扩展精度系统上才会失败。由于 log 是 Fortran 中的通用内部函数,编译器可能会全程以扩展精度格式对表达式 1.0 + x 求值,以相同的精度计算其算法,但是显然,您不能假定编译器这样做。您还可以假设一个类似示例来包含用户定义的函数。在这种情况下,编译器仍会将参数保持在扩展精度,即使函数返回了单精度结果也是如此,不过极少数(如果有)现有 Fortran 编译器也会这样做。因此,您可能希望通过将 1.0 + x 传递到变量来尝试确保对它进行一致的求值。不幸的是,如果您将该变量声明为 real,编译器仍可能会妨碍这样做,它会在变量第一次出现时使用以扩展精度格式保存在寄存器中的值替换,变量另一次出现时则使用存储在内存中的值替换。反之,您可能需要使用与扩展精度格式对应的类型来声明变量。标准 FORTRAN 77 未提供这样做的方法,而 Fortran 95 则提供了 SELECTED_REAL_KIND 机制来描述不同格式,它不显式要求以扩展精度格式对表达式求值的实现允许以该精度声明变量。简而言之,没有可移植的方法能够确保以标准 Fortran 编写此程序,防止会使证据失效的方法来对表达式 1.0 + x 求值。

在另一些示例中,即使当存储了每个子表达式从而以相同精度舍入时,在扩展精度系统上也会无法正常运行。原因是双重舍入。在缺省精度模式下,扩展精度系统最初将每个结果舍入到扩展双精度。如果随后以双精度存储该结果,则再次舍入。这两种舍入的组合生成的值与将第一个结果正确舍入到双精度获得的结果不同。如果将结果舍入到扩展双精度是“中间情况”,即正好处于两个双精度数的中间,则可能出现这种情况,因此第二次舍入由“舍入到偶数优先”规则确定。如果第二次舍入与第一次舍入方向相同,则净舍入误差将超过最后位置单元的一半。不过请注意,双重舍入只影响双精度计算。您可以证明,两个 p 位数的和、差、积或商,或者一个 p 位数的平方根在首先舍入到 q 位然后舍入到 p 位时,只要 q ≥ 2p + 2,就能得到相同的值,就像只舍入一次到 p 位一样。因此,扩展双精度足够宽,单精度计算不会经过双重舍入。

一些依赖于正确舍入的算法在双重舍入下可能会失败。实际上,即使一些算法不需要正确舍入,并且能够在各种不符合 IEEE 754 标准的计算机上正确工作,使用双重舍入时也会失败。这些方法中最有用的是可移植算法,用于执行 Goldberg 论文定理 5 中提到的模拟多精度运算。例如,定理 6 中所述用于将浮点数拆分为高位和低位的过程在双重舍入运算中无法正确运行:尝试将双精度数字 252 + 3 × 226 – 1 拆分为两部分,每一部分最多 26 位。每次运算正确舍入到双精度时,高位部分为 252 + 227,低位部分为 226 – 1,但每次运算首次舍入到扩展双精度然后再舍入到双精度时,该过程生成的高位为 252 + 228,低位为 –226 – 1。后一个数字将占据 27 位,因此在双精度下无法准确计算其平方。在扩展双精度下仍可能计算此数字的平方,但是得到的算法无法再移植到单/双精度系统上。此外,在多精度乘法算法的后面步骤中,假设所有部分乘积已经以双精度格式计算。正确处理混合的双精度和扩展双精度变量会大幅增加实现的成本。

与此类似,用于添加多精度数字(以双精度数字数组形式表示)的可移植算法也可以归类到双重舍入运算中。这些算法通常依赖类似于 Kahan 求和公式的技术。在 Goldberg 论文的 "Errors In Summation" 一节中给出的求和公式非正式说明建议,如果 sy 为浮点变量,且 |s| ≥ |y|,则可以计算如下:

t = s + y;
e = (s - t) + y;

那么,在大多数运算中,e 可以精确恢复在计算 t 时的舍入误差。此技术不适用于双重舍入运算,但是,如果 s = 252 + 1 并且 y = 1/2 – 2-54,则 s + y 以扩展双精度格式将前者舍入到 252 + 3/2,此值以双精度格式按照“舍入到偶数优先”规则舍入为 252 + 2;因此计算 t 时的净舍入误差为 1/2 + 2-54,无法准确以双精度表示,因此以上显示的表达式无法准确计算。这里再次说明,通过以扩展双精度格式计算总和,有可能弥补舍入误差,但是,程序随后需要完成额外操作以将最终输出约简回双精度,而双重舍入也会扰乱此过程。因此,虽然采用这些方法模拟多重精度运算的可移植程序可在广泛的计算机上高效正确地运行,但在扩展精度系统上则不能如宣称的这样运行。

最后,直觉上依赖于正确舍入的一些算法事实上可能会在双重舍入下正确运行。在这些情况下,处理双重舍入的成本不在于实现,而在于验证算法按照宣称方式工作。为了说明这一点,现在证明定理 7 的以下变体:

D.1.2.1 定理 7

如果 mn 是 IEEE 754 双精度模式下可表示的整数,且 |m| < 252n 具有特殊格式 n = 2i + 2j,则 (mn) ⊗ n = m,前提是假设所有浮点运算是正确舍入到双精度或者首先舍入到扩展双精度,然后舍入到双精度。


注 -  在此定理中,⊘ 和 ⊗ 分别表示计算的除法和计算的乘法。

D.1.2.2 证明

假定在无损的情况下 m > 0。让 q = mn。按二的幂换算,我们可以考虑等同的设置,其中 252m < 253,对 q 也类似,因此 mq 均为整数,最低有效位占据了单元位置(即,ulp(m) = ulp(q) = 1)。在换算之前,我们假设 m < 252,因此在扩展后,m 是偶数。此外,由于 mq 的换算值满足 m/2 < q < 2mn 的对应值必须具有两种形式之一,具体取决于 mq 哪个更大:如果 q < m,则很明显 1 < n < 2,并且由于 n 是两个二的幂之和,对于一些 kn = 1 + 2-k;与此类似,如果 q > m,则 1/2 < n < 1,因此 n = 1/2 + 2-(k+1)。由于 n 是两个二的幂之和,n 最接近于一的可能值是 n = 1 + 2-52。由于 m/(1 + 2-52) 不比小于 m 的下一个较小双精度数大,我们不能让 q = m

e 表示计算 q 时的舍入误差,因此 q = m/n + e,计算值 qn 将为 m + ne 的(一次或两次)舍入值。请首先考虑每种浮点运算正确舍入到双精度的情况。在这种情况下,|e| < 1/2。如果 n 的形式为 1/2 + 2-(k+1),则 ne = nqm 为 2-(k+1) 的整数倍,并且 |ne| < 1/4 + 2-(k+1)。这表示 |ne| ≤ 1/4。请记住,m 与下一个可表示较大数之间的差为 1,m 与下一个可表示较小数之间的差在 m > 252 时为 1,在 m = 252 时为 1/2。因此,由于 |ne| ≤ 1/4,m + ne 将舍入到 m。(即使 m = 252ne = –1/4,积仍然按照“舍入到偶数优先”规则舍入为 m。)与此类似,如果 n 的形式为 1 + 2-k,则 ne 为 2-k 的整数倍,并且 |ne| < 1/2 + 2-(k+1);这表示 |ne| ≤ 1/2。这种情况下不能让 m = 252,因为 m 严格大于 q,因此 m 与其最接近的可表示相邻数相差 ∓1。因此,由于 |ne| ≤ 1/2,m + ne 将再次舍入到 m。(即使 |ne| = 1/2,积仍然将按照“舍入到偶数优先”规则舍入为 m,因为 m 是偶数。)这就完成了正确舍入运算的证明。

在双重舍入运算中,可能仍会出现 q 是正确舍入的商(即使它实际上舍入了两次),因此如上所述,|e| < 1/2。在这种情况下,考虑到 qn 将舍入两次这样一个事实,可以采用上一段中提供的参数。为了说明这一情况,请注意,IEEE 标准要求扩展双精度格式进位至少为 64 个有效位,这样数字 m ∓ 1/2 和 m ∓ 1/4 在扩展双精度下可准确表示。因此,如果 n 的形式为 1/2 + 2-(k+1),这样 |ne| ≤ 1/4,则将 m + ne 舍入到扩展双精度一定会得到与 m 最多相差 1/4 的结果,如上所述,此值将以双精度舍入到 m。与此类似,如果 n 的形式为 1 + 2-k,这样 |ne| ≤ 1/2,则将 m + ne 舍入到扩展双精度一定会得到与 m 最多相差 1/2 的结果,此值将以双精度舍入到 m。请记住,这种情况中的 m > 252

最后,我们只剩下考虑由于双重舍入而使 q 不是正确的舍入商的情况。在这些情况下,最坏的情况下将得到 |e| < 1/2 + 2-(d + 1),其中 d 是扩展双精度中额外位的数量。所有现有的扩展精度系统支持恰好具有 64 个有效位的扩展双精度格式;对于此格式,d = 64 – 53 = 11。如果第二次舍入由“舍入到偶数优先”规则确定,则双重舍入只会得到不正确的舍入结果,因而 q 必须为偶数。因此,如果 n 的形式为 1/2 + 2-(k + 1),则 ne = nqm 是 2-k 的整数倍,并且 |ne| < (1/2 + 2-(k + 1))(1/2 + 2-(d + 1)) = 1/4 + 2-(k + 2) + 2-(d + 2) + 2-(k + d + 2)。

如果 kd,则表明 |ne| ≤ 1/4。如果 k > d,就有了 |ne| ≤ 1/4 + 2-(d + 2)。在任何一种情况下,积的第一次舍入将得到与 m 最多相差 1/4 的结果,按照以前的参数,第二次舍入将舍入到 m。与此类似,如果 n 的形式为 1 + 2-k,则 ne 为 2-(k - 1) 的整数倍,并且 |ne| < 1/2 + 2-(k + 1) + 2-(d + 1) + 2-(k + d + 1)

如果 kd,则表明 |ne| ≤ 1/2。如果 k > d,就有了 |ne| ≤ 1/4 + 2-(d + 1)。在任何一种情况下,积的第一次舍入将得到与 m 最多相差 1/2 的结果,仍旧按照以前的参数,第二次舍入将舍入到 m。 ♦

以上证据表明,只有在商进行双重舍入时,积才会发生双重舍入。即使这样,也会舍入到正确结果。该证据还表明,扩展我们的推理以包括双重舍入可能会遇到难题,即使程序只有两次浮点运算也是如此。对于更复杂的程序,可能无法系统地解释双重舍入效果,更不用说那些更常见的双精度和扩展双精度计算组合。

D.1.3 扩展精度的编程语言支持

以上示例并非用于说明扩展精度本身有害。当程序员具备选择性使用的能力时,许多程序可从扩展精度中获益。不幸的是,当前编程语言没有提供足够的手段,来方便程序员指定适合使用扩展精度的情况和方法。为了指明需要什么支持,请考虑您希望用于管理扩展精度使用的方式。

在使用双精度作为其名义工作精度的可移植程序中,可以采用五种方法来控制更宽精度的使用:

  1. 在扩展精度系统上尽可能使用扩展精度进行编译,这样可生成最快的代码。很明显,大多数数值软件所需要的算术精度并不会超过每次运算中由“计算机精度”所限制的相对误差。当内存中的数据以双精度存储时,计算机精度通常采用该双精度的最大相对舍入误差,因为假定在输入数据时这些数据已进行舍入,并且运算结果将在存储时同样地进行舍入。因此,虽然以扩展精度计算某些中间结果可能会产生更准确的结果,但扩展精度并不是必要的。在这种情况下,我们可能更希望仅当编译器使用扩展精度不会明显减慢程序速度时,才使用扩展精度,否则使用双精度。

  2. 如果比双精度更宽的某种格式明显更快并且足够宽,则使用该格式,否则采用其他格式。当扩展精度可用时,可以更轻松地进行某些计算,不过这些计算也可以采用双精度运行,只需稍微付出一些劳动。请考虑计算双精度数字的欧几里得向量范数。通过使用 IEEE 754 扩展双精度格式及其更宽的指数范围来计算元素的平方并累积计算其总和,我们可以轻松地避免可行长度的向量过早下溢或溢出。在扩展精度系统上,这是计算范数的最快方法。在单/双精度系统上,如果扩展双精度格式根本不受支持,则必须在软件中进行模拟,并且此类模拟将比仅使用双精度要慢得多。需要测试异常标志以确定是否发生下溢或上溢,如果发生,将使用显式缩放重复计算过程。请注意,为支持扩展精度的这种用法,一种语言必须指示相当快的最宽可用格式(以便程序可以选择使用哪个方法),还必须提供环境参数来指示每种格式的精度和范围(以便程序可以验证相当快的最宽格式足够宽,例如,它比双精度的范围更宽)。

  3. 使用比双精度更宽的某种格式,即使该格式必须在软件中进行模拟。对于比欧几里得范数示例更复杂的程序,程序员可能希望避免编写两种版本的程序,转而依赖于扩展精度(即使它速度比较慢)。另外,语言必须提供环境参数,这样程序可以确定最宽可用格式的范围和精度。

  4. 不使用更宽的精度;将结果正确地舍入到双精度格式的精度,即使它可能具有更宽的范围。最容易编写的程序取决于正确舍入的双精度运算,包括上面提到的一些示例。语言必须为程序员提供一种方式来指示不得使用扩展精度,即使可在寄存器中使用比双精度更宽的指数范围来计算中间结果也不行。对于使用这种方式计算的中间结果,如果在存储到内存时发生下溢,仍然可能执行双重舍入:如果算术运算结果首先舍入为 53 个有效位,然后在必须非规范化时再舍入为更少的有效位,最终结果可能不同于仅舍入一次为非规范化数字所得到的结果。当然,这种形式的双重舍入几乎不可能对任何实际程序产生负面影响。

  5. 将结果正确地舍入为双精度格式的精度和范围。这种双精度的严格实现最有益于下面这样的程序:这些程序在双精度格式的精度和范围限制附近测试数字软件或算术本身。这类严谨的测试程序一般很难以可移植方式编写;当它们必须利用伪子例程和其他技巧强制将结果舍入为特定格式时,它们甚至更难于编写并且容易出错。因此,程序员在使用扩展精度系统来开发必须可移植到所有 IEEE 754 实现的强大软件时,很快会欣喜地发现,能够模仿单/双精度系统的运算且不需要付出额外工作。

当前没有语言全部支持上面所有五种选项。实际上,少量语言已经尝试使程序员能够从根本上控制对扩展精度的使用。一个值得注意的例外情况是 ISO/IEC 9899:1999 编程语言 - C 标准,它是对 C 语言的主修订版。

C99 标准允许实现在表达式求值时使用的格式比通常与表达式类型关联的格式更宽,但 C99 标准推荐使用仅有的三种表达式求值方法中的一种。这三种推荐方法的特征在一定程度上为哪些表达式“被提升”为更宽的格式。通过定义预处理程序宏 FLT_EVAL_METHOD 来帮助实现方式确定使用的是哪种方法:如果 FLT_EVAL_METHOD 为 0,则以每个表达式的类型所对应的格式对表达式求值;如果 FLT_EVAL_METHOD 为 1,则将 float 表达式提升为对应于 double 的格式;如果 FLT_EVAL_METHOD 为 2,则将 floatdouble 表达式提升为对应于 long double 的格式。(允许实现将 FLT_EVAL_METHOD 设置为 –1,以指示表达式求值方法是不确定的。)C99 标准还要求在 <math.h> 头文件中定义类型 float_tdouble_t。这两种类型相应地至少与 floatdouble 一样宽,其目的是为了与 floatdouble 表达式的求值类型相符。例如,如果 FLT_EVAL_METHOD 为 2,则 float_tdouble_tlong double。最后,C99 标准要求在 <float.h> 头文件中定义预处理程序宏,以指定对应于每种浮点类型的格式的范围和精度。

C99 标准要求或推荐的功能组合支持上面所列五个选项中的一部分,但并非全部支持。例如,如果实现将 long double 类型映射到扩展双精度格式,并将 FLT_EVAL_METHOD 定义为 2,则程序员可以合理地假定扩展精度相对较快,因此欧几里得范数示例这样的程序只需使用 long double 类型(或 double_t)的中间变量。另一方面,即使匿名表达式存储在内存中(例如,当编译器必须使浮点寄存器溢出时),同一实现也必须将匿名表达式保留为扩展精度,并且必须存储赋值给变量(这些变量已声明 double 来转换为双精度)的表达式结果,即使这些结果可能已保留在寄存器中。因此,doubledouble_t 类型都不能编译为在当前扩展精度硬件上生成最快的代码。

类似地,对于本节中的示例所阐述的一些问题(但不是全部),C99 标准提供了解决方法。如果将表达式 1.0 + x 赋值给变量(任何类型)并且到处使用该变量,log1p 函数的 C99 标准版本保证可以正确工作。但是,编写可移植的高效 C99 标准程序来将双精度数字拆分为高位部分和低位部分会更加困难:如果我们无法保证 double 表达式正确舍入为双精度,如何在正确位置拆分并避免双重舍入?一个解决方法是使用 double_t 类型在单/双精度系统上拆分双精度,并在扩展精度系统上拆分扩展精度,以便在其中任一情况下都能正确舍入运算。在定理 14 中说明,假设我们知道基础运算的精度,则我们可以在任何位的位置进行拆分;并且 FLT_EVAL_METHOD 和环境参数宏应该向我们提供此信息。

下面的代码片段显示了一个可能的实现:

#include <math.h>
#include <float.h>
 
#if (FLT_EVAL_METHOD==2)
#define PWR2  LDBL_MANT_DIG - (DBL_MANT_DIG/2)
#elif ((FLT_EVAL_METHOD==1) || (FLT_EVAL_METHOD==0))
#define PWR2  DBL_MANT_DIG - (DBL_MANT_DIG/2)
#else
#error FLT_EVAL_METHOD unknown!
#endif
 
...
    double   x, xh, xl;
    double_t m;
 
    m = scalbn(1.0, PWR2) + 1.0;  // 2**PWR2 + 1
    xh = (m * x) - ((m * x) - x);
    xl = x - xh;

要找到此解决方法,您必须了解:可以采用扩展精度对 double 表达式求值,随后发生的双重舍入问题会导致算法出现问题,并且可以根据定理 14 改用扩展精度。一个更明显的解决方法是,只需指定将每个表达式正确舍入为双精度。在扩展精度系统上,这样做只需要更改舍入精度模式,但很不幸,C99 标准并未提供可移植方法来完成此操作。在浮点 C 修订的早期草案工作文档中,指定了为支持浮点而对 C90 标准进行的更改,并建议系统上具有舍入精度模式的实现提供 fegetprecfesetprec 函数来获取和设置舍入精度,这类似于获取和设置舍入方向的 fegetroundfesetround 函数。对 C99 标准进行更改之前,已经删除了此建议。

巧合的是,在 C99 标准用于支持在具有不同整数运算功能的系统之间进行移植的方法中,给出了一种更好的方法来支持不同的浮点体系结构。每个 C99 标准实现提供一个 <stdint.h> 头文件,其中定义该实现支持的那些整数类型,这些类型根据其大小和效率进行命名:例如,int32_t 是一种恰好为 32 位宽度的整数类型,int_fast16_t 是该实现的最快整数类型(至少为 16 位宽度),而 intmax_t 是能够支持的最宽整数类型。可以想像对于浮点类型有一个类似的方案:例如,float53_t 可以命名恰好为 53 位精度但范围可能更宽的浮点类型,float_fast24_t 可以命名实现中至少为 24 位精度的最快类型,而 floatmax_t 可以命名能够支持的最宽且相当快的类型。在扩展精度系统上,快类型可以允许编译器生成可能最快的代码,只是命名变量值不得在寄存器溢出后更改。在扩展精度系统上,精确宽度类型将导致编译器将舍入精度模式设置为舍入到指定精度,从而允许在相同约束限制下的更宽范围。最后,double_t 可以使用 IEEE 754 双精度格式的精度和范围来命名类型,从而提供严格的双精度求值。随着相应地命名环境参数宏,此类方案将轻松支持上面描述的所有五种选项,并允许程序员轻松地明确指示其程序需要的浮点语义。

对扩展精度的语言支持一定如此复杂吗?在单/双精度系统上,上面所列五种选项中的四种都一致,并且不需要区分快类型和精确宽度类型。但是,扩展精度系统带来了困难的选择:此类系统支持的纯双精度计算或纯扩展精度计算的效率不如两者的混合形式,而不同程序要求不同的混合形式。此外,使用扩展精度的选择时机不应留给编译器编写人员,他们经常受到基准的迷惑,有时数字分析员会直率地告诉他们浮点运算被视为“天生不精确”,因此整数运算既不值得预测,也无法预测。相反,选择权必须交给程序员,他们需要使用能够表达其选择的语言。

D.1.4 小结

上述评论并不是为了贬低扩展精度系统,而是揭示几个谬误。第一个谬误是所有 IEEE 754 系统针对同一程序必须产生相同的结果。我们已经着重讨论了扩展精度系统与单/双精度系统之间的区别,但在这些系统之间,每个系列仍然有进一步的差别。例如,一些单/双精度系统提供了一个指令,将两个数字相乘,然后加第三个数字,最后只进行一次舍入。此运算称为混合乘加运算,可导致同一程序在不同单/双精度系统上生成不同的结果。另外,类似于扩展精度,它甚至会导致同一程序在同一系统上生成不同结果,具体取决于是否使用此运算以及使用时间。虽然可以通过不可移植方式使用混合乘加运算来执行多精度乘法而无需拆分,但该运算还会阻挠定理 6 的拆分过程。即使 IEEE 标准未预见到此类运算,但它仍然符合要求:中间乘积提交给不受用户控制的“目标”,该目标的宽度足以精确地保存中间结果,最终的总和将正确地舍入,以符合其单精度或双精度目标。

IEEE 754 精确地规定指定程序必须提交的结果的理念仍然是非常吸引人的。许多程序员认为,他们了解程序的行为,并且可以证明程序将正确工作,而不用考虑编译该程序的编译器或运行该程序的计算机。在许多情况下,支持这一信念正是计算机系统和编程语言的设计者努力的目标。不幸的是,当遇到浮点运算时,这一目标实际上无法实现。IEEE 标准的作者知道这一点,并且没有尝试实现这一目标。因此,在整个计算机行业中,尽管几乎普遍遵守大多数 IEEE 754 标准,但可移植软件的程序员必须继续应付不可预测的浮点运算。

如果程序员打算利用 IEEE 754 的功能,他们需要使用能够使浮点运算变得可预测的编程语言。C99 标准在一定程度上改善了可预测性,代价是需要程序员编写程序的多个版本,对每个 FLT_EVAL_METHOD 编写一个版本。在未来的语言中,是否改为选择允许程序员编写单个程序,并且使用能够明确表示对 IEEE 754 语义的依赖程度的语法,这一点仍可期待。从现有扩展精度系统上,我们无法预见到这种远景,因为它们促使我们假定编译器和硬件比程序员更了解在给定系统上执行计算的方式。该假定是第二个谬误:计算结果中需要的精确度不依赖于产生该结果的计算机,而只依赖于从该结果得到的结论。对于程序员、编译器和硬件,最多只有程序员可能了解这些结论是什么。