oracle解一元二次方程,第 6 章浮点运算-CSDN博客

第 6 章浮点运算

本章分析浮点运算并提出避免和检测数值计算错误的策略。

有关 SPARC 和 x86 处理器浮点计算的详细说明，参见《数值计算指南》。

6.1 简介

SPARC 处理器上的 Fortran 95 浮点环境实现了“二进制浮点运算 IEEE 标准 754”指定的运算模型。该环境使您能够开发强健、高性能、可移植的数值应用程序。它还提供了用于分析研究数值程序任何不正常行为的工具。

在数值程序中，有许多潜在因素可引起计算错误：

计算模型错误

所使用的算法在数值上不稳定

病态数据

硬件产生意外结果

找出数值计算中出错的原因非常困难。可以尽可能使用市面上销售的经过测试的库程序包来减少编码错误几率。算法的选择是另一个关键问题。使用合适的计算机运算同样也是一个关键问题。

本章不打算教授或解释数值错误分析。此处提供的资料旨在介绍 Fortran 95 实现的 IEEE 浮点模型。

6.2 IEEE 浮点运算

IEEE 运算是一种处理会导致如下问题的运算操作的相对较新的方法：无效操作数、被零除、上溢、下溢或结果不精确。不同之处在于舍入、近零数字的处理以及接近机器最大值的数字的处理。

IEEE 标准支持用户处理异常、舍入和精度。因此，此标准支持区间运算和异常诊断。IEEE 标准 754 可以使诸如 exp 和 cos 之类的基本函数标准化，以便创建高精确度的运算，并结合数值和符号代数计算。

与其他任何种类的浮点运算相比，IEEE 运算向用户提供了对计算的更强大控制。此标准简化了编写复杂可移植数值程序的任务。很多有关浮点算法的问题都涉及数的基本运算。例如：

当计算机硬件无法表示无限精确的结果时，运算结果会怎样？

有些基本运算(比如乘法和加法)是否具有交换性？

另一类问题与浮点异常及异常处理有关。如果进行下列运算会发生什么情况：

二个非常大的同号数字相乘？

非零值除以零？

零除以零？

在较早的运算模型中，第一类问题可能不会得到预期的答案，而第二类问题中的异常情况可能都会得到相同的结果：程序立即中止或使用无用结果继续执行。

此标准可确保运算产生具有预期性质的、符合数学预期的结果。它还确保异常情况产生指定的结果，除非用户明确地指定其他选项。

例如，直观地引入了异常值 +Inf、-Inf 和 NaN：

big*big = +Inf 正无穷

big*(-big) = -Inf 负无穷

num/0.0 = +Inf 其中 num

> 0.0

num/0.0 = -Inf 其中 num

< 0.0

0.0/0.0 = NaN 非数字

并且，还标识出五种类型的浮点异常：

无效。运算操作数在数学上无效－例如，0.0/0.0、sqrt(-1.0) 和 log(-37.8)

被零除。除数为零，被除数为有限的非零数字－例如，9.9/0.0

上溢。运算产生的结果超出指数范围－例如，MAXDOUBLE+0.0000000000001e308

下溢。运算产生的结果太小，无法用正常数字表示－例如，MINDOUBLE * MINDOUBLE

不精确。运算产生的结果无法用无限精度来表示－例如，输入中的 2.0/3.0、log(1.1) 和 0.1

在《数值计算指南》中介绍了 IEEE 标准的实现。

6.2.1 –ftrap=mode 编译器选项

-ftrap=mode 选项可捕获浮点异常。如果 ieee_handler() 调用未建立信号处理程序，异常会用内存转储核心转储文件终止程序。有关该编译器选项的详细信息，请参见《Fortran 用户指南》。例如，为了能够捕获上溢、被零除和无效运算，可使用 -ftrap=common 进行编译。(这是 f95 的缺省选项。)

注 –

必须使用 -ftrap= 编译应用程序的主程序才能进行捕获。

6.2.2 浮点异常

f95 程序不会自动报告异常。要显示程序终止时产生的浮点异常列表，需要显式调用 ieee_retrospective(3M)。一般情况下，如果发生了无效、被零除或上溢异常中的任何一种，都会产生消息。不精确异常不会产生消息，因为它们在实际程序中发生得过于频繁。

6.2.2.1 回顾性摘要

ieee_retrospective函数对浮点状态寄存器进行查询，以找出产生了哪些异常，并打印一条有关标准错误的消息来通知您哪些是已产生但尚未清除的异常。此消息通常如下所示；格式可能会随各编译器版本而变化：

Note: IEEE floating-point exception flags raised:

Division by Zero;

IEEE floating-point exception traps enabled:

inexact; underflow; overflow; invalid operation;

See the Numerical Computation Guide, ieee_flags(3M),

ieee_handler(3M)

Fortran 95 程序需要显式调用 ieee_retrospective，并使用 -xlang=f77 进行编译，以便与 f77 兼容库进行链接。

用 -f77 兼容标志进行编译，将启用程序终止时自动调用 ieee_retrospective 惯例。

可以使用 ieee_flags() 关闭任意或所有这些消息，方法是在调用 ieee_retrospective 之前清除异常状态标记。

6.2.3 处理异常

根据 IEEE 标准进行异常处理是 SPARC 和 x86 处理器上的缺省异常处理方法。但是，检测浮点异常和生成浮点异常信号(SIGFPE)之间是有区别的。

按照 IEEE 标准，当浮点运算期间出现未捕获的异常时，会发生两件事情：

系统返回缺省结果。例如，对于 0/0 (invalid)，系统返回结果为 NaN。

会设置标志来指示引起了异常。例如，对于 0/0 (invalid)，系统会设置“无效运算”标志。

6.2.4 捕获浮点异常

f95 在处理浮点异常方面与早期的 f77 编译器有着明显的区别。

缺省情况下，f95 会自动捕获被零除、上溢和无效运算。而 f77 在缺省情况下不会为浮点异常自动产生信号来中断正在运行的程序。其假设是：只要返回期望的值，大多数异常都无关紧要，对其进行捕获会降低性能。

可以使用 f95 命令行选项 -ftrap 来更改缺省设置。f95 的缺省选项是 -ftrap=common。要按早期的 f77 缺省设置进行运算，请使用 -ftrap=%none 编译主程序。

6.2.5 非标准运算

可以手动禁用标准 IEEE 运算中一个称为渐进下溢的特征。禁用后，程序将被视为在非标准运算下运行。

IEEE 运算标准规定了一种通过动态调整有效数的小数点来渐进处理下溢结果的方法。按 IEEE 浮点格式，小数点出现在有效数之前，并且有一个隐式前导位 1。当浮点计算结果会下溢时，渐进下溢允许将此隐式前导位清为 0，并将小数点移入有效数之中，否则，浮点计算结果便会产生下溢。对于 SPARC 处理器，该结果不是在硬件而是在软件中完成的。如果程序产生的下溢很多(也许这表示您的算法有问题)，可能会导致性能损失。

可以通过以下方式禁用然后关闭渐进下溢：使用 -fns 选项进行编译，或从程序内调用库例程 nonstandard_arithmetic()。调用 standard_arithmetic() 可重新开启渐进下溢。

注 –

为提高效率，必须用 -fns 编译应用程序的主程序。参见《Fortran 用户指南》。

对于传统应用程序，请注意：

standard_arithmetic() 子例程会替换名为 gradual_underflow() 的早期例程。

nonstandard_arithmetic() 子例程会替换一个早期例程(名为 abrupt_underflow())。

注 –

-fns 选项和 nonstandard_arithmetic() 库例程仅在某些 SPARC 系统中有效。

6.3 IEEE 例程

以下接口可帮助用户使用 IEEE 运算，并在手册页中进行说明。这些接口多数都在数学库 libsunmath 和几个 .h 文件中。

ieee_flags(3m)－控制舍入方向和舍入精度；查询异常状态；清除异常状态

ieee_handler(3m)－建立异常处理例程

ieee_functions(3m)－列出每个 IEEE 函数的名称和用途

ieee_values(3m)－列出返回特殊值的函数

本部分介绍的其他 libm 函数：

ieee_retrospective

nonstandard_arithmetic

standard_arithmetic

SPARC 处理器对不同方面的软硬件支持组合符合 IEEE 标准。

最新的 SPARC 处理器包含具有整数乘法和除法指令以及硬件平方根的浮点单元。

当编译代码与运行时浮点硬件正确匹配时，会获得最佳性能。编译器的 -xtarget= 选项允许指明运行时硬件。例如，-xtarget=ultra 会通知编译器生成在 UltraSPARC 处理器上执行效果最佳的目标代码。

实用程序 fpversion 显示安装了哪种浮点硬件，并指示要指定的合适的 -xtarget 值。该实用程序可在所有 Sun SPARC 体系结构中运行。有关详细信息，请参见 fpversion(1)、《Fortran 用户指南》和《数值计算指南》。

6.3.1 标志和 ieee_flags()

ieee_flags 函数用于查询和清除异常状态标志。它是 Sun 编译器随带的 libsunmath 库的一部分，可执行下列任务：

控制舍入方向和舍入精度

检查异常标志状态

清除异常状态标志

ieee_flags 调用的一般形式为：

flags = ieee_flags( action, mode, in, out )

四个参数中的每一个都是字符串。输入为 action、mode 和 in。输出为 out 和 flags。ieee_flags 是一个整数值函数。flags 中返回有用的信息，作为 1 位标志集合。有关完整信息，请参见 ieee_flags(3m) 手册页。

下表显示了所有可能的参数值。

表 6–1 ieee_flags ( action, mode, in, out ) 参数值参数允许值

actionget, set, clear, clearall

modedirection, exception

in, outnearest, tozero, negative, positive, extended, double single, inexact, division, underflow, overflow, invalid all, common

注意，这些是文字字符串，且输出参数 out 必须至少是 CHARACTER*9。in 和 out 的可能值的含义取决于与其一起使用的 action 和 mode。下表对此进行了概括：

表 6–2 ieee_flags in、 out 参数的含义in 和 out 的值所指

nearest, tozero, negative, positive舍入方向

extended, double, single舍入精度

inexact, division, underflow, overflow, invalid异常

all全部五种异常

common常见异常：无效、除法、上溢

例如，要确定引起了标志的具有最高优先级的异常，请将输入参数 in 作为空字符串传递：

CHARACTER *9, out

ieeer = ieee_flags( ’get’, ’exception’, ’’, out )

PRINT *, out, ’ flag raised’

另外，要确定是否引起了 overflow 异常标志，请将输入参数 in 设置为 overflow。返回时，如果 out 等于 overflow，会出现 overflow 异常标志；否则不会出现该标志。

ieeer = ieee_flags( ’get’, ’exception’, ’overflow’, out )

IF ( out.eq. ’overflow’) PRINT *,’overflow flag raised’

示例：清除 invalid 异常：

ieeer = ieee_flags( ’clear’, ’exception’, ’invalid’, out )

示例：清除所有异常：

ieeer = ieee_flags( ’clear’, ’exception’, ’all’, out )

示例：将舍入方向设置为零：

ieeer = ieee_flags( ’set’, ’direction’, ’tozero’, out )

示例：将舍入精度设置为 double：

ieeer = ieee_flags( ’set’, ’precision’, ’double’, out )

6.3.1.1 用 ieee_flags 关闭所有警告消息。

使用清除 action 调用 ieee_flags(如下例所示)可以重置任何未清除的异常。在程序退出之前进行该调用，可禁止系统在程序终止时产生浮点异常警告消息。

示例：用 ieee_flags() 清除所有产生的异常：

i = ieee_flags(’clear’, ’exception’, ’all’, out )

6.3.1.2 用 ieee_flags 检测异常

以下示例演示如何确定早期计算引起的浮点异常。会将系统 include 文件 floatingpoint.h 中定义的位屏蔽应用于 ieee_flags 的返回值。

在以下示例(即 DetExcFlg.F)中，include 文件是使用 #include 预处理程序指令引入的，这就要求以 .F 后缀命名源文件。下溢是由最小的双精度数除以 2 引起的。

示例：使用 ieee_flags 检测异常，然后对其进行解码：

#include "floatingpoint.h"

CHARACTER*16 out

DOUBLE PRECISION d_max_subnormal, x

INTEGER div, flgs, inv, inx, over, under

x = d_max_subnormal() / 2.0 ! Cause underflow

flgs=ieee_flags(’get’,’exception’,’’,out) ! Which are raised?

inx = and(rshift(flgs, fp_inexact) , 1) ! Decode

div = and(rshift(flgs, fp_division) , 1) ! the value

under = and(rshift(flgs, fp_underflow), 1) ! returned

over = and(rshift(flgs, fp_overflow) , 1) ! by

inv = and(rshift(flgs, fp_invalid) , 1) ! ieee_flags

PRINT *, "Highest priority exception is: ", out

PRINT *, ’ invalid divide overflo underflo inexact’

PRINT ’(5i8)’, inv, div, over, under, inx

PRINT *, ’(1 = exception is raised; 0 = it is not)’

i = ieee_flags(’clear’, ’exception’, ’all’, out) ! Clear all

END

示例：编译并运行上述示例 (DetExcFlg.F)：

demo% f95 DetExcFlg.F

demo% a.out

Highest priority exception is: underflow

invalid divide overflo underflo inexact

0 0 0 1 1

(1 = exception is raised; 0 = it is not)

demo%

6.3.2 IEEE 极值函数

编译器提供了一个函数集，可以调用其中的函数来返回特殊的 IEEE 极值。这些值，如 infinity 或 minimum normal，可以直接在应用程序中使用。

示例：基于硬件支持的最小数值的收敛测试如下所示：

IF ( delta .LE. r_min_normal() ) RETURN

下表列出了可用的值：

表 6–3 返回 IEEE 值的函数IEEE 值双精度单精度

infinityd_infinity()r_infinity()

quiet NaNd_quiet_nan()r_quiet_nan()

signaling NaNd_signaling_nan()r_signaling_nan()

min normald_min_normal()r_min_normal()

min subnormald_min_subnormal()r_min_subnormal()

max subnormald_max_subnormal()r_max_subnormal()

max normald_max_normal()r_max_normal()

两个 NaN 值(quiet 和 signaling)是无序的，不能用于比较，如 IF(X.ne.r_quiet_nan())THEN...。要确定某些值是否是 NaN，请使用函数 ir_isnan(r) 或 id_isnan(d)。

以下手册页列出了这些函数的 Fortran 名称：

libm_double(3f)

libm_single(3f)

ieee_functions(3m)

另请参见：

ieee_values(3m)

floatingpoint.h 头文件和 floatingpoint(3f)

6.3.3 异常处理程序和 ieee_handler()

关于 IEEE 异常，通常需要关注以下问题：

异常出现时会发生什么情况？

如何使用 ieee_handler() 来建立一个可用作异常处理程序的用户函数？

如何编写可用作异常处理程序的函数？

如何定位异常－即异常出现在何处？

用户例程的异常捕获以系统产生浮点异常信号开始。signal:

floating-point exception的 UNIX 标准名称是 SIGFPE。出现异常时，SPARC 平台上的缺省情况是不产生 SIGFPE。要使系统产生 SIGFPE，必须先启用异常捕获，这通常通过对 ieee_handler() 的调用来完成。

6.3.3.1 建立异常处理程序函数

要将函数作为异常处理程序建立，请将函数名称与要监视的异常的名称和要采取的操作一起传递给 ieee_handler()。一旦建立了处理程序，无论何时出现特定的浮点异常和调用指定的函数，都会产生 SIGFPE 信号。

ieee_handler() 的调用形式如下表所示：

表 6–4 ieee_handler (action , exception , handler) 的参数参数类型可能值

actioncharacterget、set 或 clear

exceptioncharacterinvalid、division、overflow、underflow 或 inexact

handler函数名用户处理函数的名称或 SIGFPE_DEFAULT、SIGFPE_IGNORE 或 SIGFPE_ABORT

返回值integer0 =OK

用 f95 编译的、调用 ieee_handler() 的 Fortran 95 例程还应该声明：

#include ’floatingpoint.h’

特殊参数 SIGFPE_DEFAULT、SIGFPE_IGNORE 和 SIGFPE_ABORT 定义在这些包含文件中，可用于更改与特定异常相应的程序行为：

SIGFPE_DEFAULT 或 SIGFPE_IGNORE出现指定异常时不采取任何操作。

SIGFPE_ABORT程序在异常时中止(可能会使用转储文件)。

6.3.3.2 编写用户异常处理程序函数

异常处理程序采取的操作由您决定。但是，此例程必须是整型函数，且具有下面指定的三个参数：

handler_name( sig, sip, uap )

handler_name 是此整型函数的名称。

sig 是一个整数。

sip 是具有结构 siginfo 的记录。

未使用 uap。

示例：异常处理程序函数：

INTEGER FUNCTION hand( sig, sip, uap )

INTEGER sig, location

STRUCTURE /fault/

INTEGER address

INTEGER trapno

END STRUCTURE

STRUCTURE /siginfo/

INTEGER si_signo

INTEGER si_code

INTEGER si_errno

RECORD /fault/ fault

END STRUCTURE

RECORD /siginfo/ sip

location = sip.fault.address

... actions you take ...

END

此示例需要修改才能运行在 64 位 SPARC 体系结构上，方法是使用 INTEGER*8 替换每个 STRUCTURE 中的所有 INTEGER 声明。

如果由 ieee_handler() 启用的处理程序例程与此示例中一样，是用 Fortran 编写的，则此例程不能对其第一个参数 (sig) 进行任何引用。该第一个参数按值传递给此例程，并且只能作为 loc(sig) 进行引用。此值是信号编号。

通过处理程序检测异常

下列示例展示如何创建处理程序例程来检测浮点异常。

示例：检测异常并中止：

demo% cat DetExcHan.f

EXTERNAL myhandler

REAL :: r = 14.2 , s = 0.0

i = ieee_handler (’set’, ’division’, myhandler )

t = r/s

END

INTEGER FUNCTION myhandler(sig,code,context)

INTEGER sig, code, context(5)

CALL abort()

END

demo% f95 DetExcHan.f

demo% a.out

Abort

demo%

SIGFPE 可在浮点异常出现的任何时间产生。检测到 SIGFPE 时，控制将传递给 myhandler 函数，该函数会立即中止。用 -g 编译，并使用 dbx 查找异常位置。

通过处理程序定位异常

示例：定位异常(打印地址)并中止：

demo% cat LocExcHan.F

#include "floatingpoint.h"

EXTERNAL Exhandler

INTEGER Exhandler, i, ieee_handler

REAL:: r = 14.2 , s = 0.0 , t

C Detect division by zero

i = ieee_handler( ’set’, ’division’, Exhandler )

t = r/s

END

INTEGER FUNCTION Exhandler( sig, sip, uap)

INTEGER sig

STRUCTURE /fault/

INTEGER address

END STRUCTURE

STRUCTURE /siginfo/

INTEGER si_signo

INTEGER si_code

INTEGER si_errno

RECORD /fault/ fault

END STRUCTURE

RECORD /siginfo/ sip

WRITE (*,10) sip.si_signo, sip.si_code, sip.fault.address

10 FORMAT(’Signal ’,i4,’ code ’,i4,’ at hex address ’, Z8 )

Exhandler=1

CALL abort()

END

demo% f95 -g LocExcHan.F

demo% a.out

Signal 8 code 3 at hex address 11230

Abort

demo%

在 64 位 SPARC 环境中，请用 INTEGER*8 替换每个 STRUCTURE 中的 INTEGER 声明，用 i8 替换 i4 格式。(注意，该例接受 VAX Fortran STRUCTURE 语句，依靠的是 f95 编译器的扩展。)

大多数情况下，知道异常的实际地址并无太大用处，但对于 dbx 除外：

demo% dbx a.out

(dbx) stopi at 0x11230 Set breakpoint at address

(2) stopi at &MAIN+0x68

(dbx) run Run program

Running: a.out

(process id 18803)

stopped in MAIN at 0x11230

MAIN+0x68: fdivs %f3, %f2, %f2

(dbx) where Shows the line number of the exception

=>[1] MAIN(), line 7 in "LocExcHan.F"

(dbx) list 7 Displays the source code line

7 t = r/s

(dbx) cont Continue after breakpoint, enter handler routine

Signal 8 code 3 at hex address 11230

abort: called

signal ABRT (Abort) in _kill at 0xef6e18a4

_kill+0x8: bgeu _kill+0x30

Current function is exhandler

24 CALL abort()

(dbx) quit

demo%

当然，还有更容易的方法来确定引起错误的源码行。但是，本例确实足以展示异常处理的基本内容。

6.4 调试 IEEE 异常

定位异常出现的位置需要启用异常捕获。这可以通过使用 -ftrap=common 选项(用 f95 编译器编译时的缺省选项)进行编译，或通过使用 ieee_handler() 建立异常处理程序例程来完成。启用了异常捕获，便可从 dbx 中运行程序，使用 dbx catch FPE 命令来查看出错位置。

使用 -ftrap=common 编译的优点是：无需修改源代码即可捕获异常。但是，通过调用 ieee_handler()，您可以更有选择性地查看异常。

示例：dbx 的编译及使用：

demo% f95 -g myprogram.f

demo% dbx a.out

Reading symbolic information for a.out

...

(dbx) catch FPE

(dbx) run

Running: a.out

(process id 19739)

signal FPE (floating point divide by zero) in MAIN at line 212 in file "myprogram.f"

212 Z = X/Y

(dbx) print Y

y = 0.0

(dbx)

如果发现程序因上溢和其他异常而终止，可调用 ieee_handler() 明确定位第一处上溢，以便只捕获上溢。这至少需要修改主程序的源代码，如下例所示。

示例：在其他异常出现时定位上溢：

demo% cat myprog.F

#include “floatingpoint.h”

program myprogram

...

ier = ieee_handler(”set’,’overflow’,SIGFPE_ABORT)

...

demo% f95 -g myprog.F

demo% dbx a.out

Reading symbolic information for a.out

...

(dbx) catch FPE

(dbx) run

Running: a.out

(process id 19793)

signal FPE (floating point overflow) in MAIN at line 55 in file "myprog.F"

55 w = rmax * 200. ! Cause of the overflow

(dbx) cont ! Continue execution to completion

execution completed, exit code is 0

(dbx)

为了具有选择性，此示例引入了 #include，它需要以 .F 后缀重命名源文件，并调用 ieee_handler()。您可更深入一层，创建出现上溢异常时要调用的自己的处理程序函数，执行一些应用程序特定的分析，并在中止前打印中间结果或调试结果。

6.5 更深层次的数值风险

本部分解决一些运算操作无意间可能会产生无效、被零除、上溢、下溢或不精确异常的实际问题。

例如，在 IEEE 标准之前，如果在计算机中将两个非常小的数相乘，结果可能为零。多数大型机和小型机的行为亦是如此。使用 IEEE 运算，渐进下溢会扩大动态计算范围。

例如，假设某一 32 位处理器以 1.0E-38 作为机器中可表示的最小值 epsilon。将两个小数相乘：

a = 1.0E-30

b = 1.0E-15

x = a * b

较早的运算会得到 0.0，但如果使用 IEEE 运算和相同的字长，却会得到 1.40130E-45。此时便出现了下溢，告诉您答案比机器自然表示的值小。该结果是通过从尾数中“窃取”一些位并将其移交给指数来完成的。得到的结果(即一个非规范化数)从某种意义上说精确度比较低，但从另外一种意义上说精确度又比较高。其深层含意已超出了本次讨论的范围。如果有兴趣，可以参考《Computer》1980 年 1 月，第 13 卷，第 1 期，尤其是 J. Coonen 的文章 "Underflow

and the Denormalized Numbers"。

大多数科学程序都有对舍入很敏感的代码段，通常是在方程求解或矩阵因子分解中。若不采用渐进下溢，程序员就得自己实现检测接近不准确阈值的方法。否则，他们就必须放弃对实现强大、稳定算法的追求。

有关这些主题的详细信息，请参见《数值计算指南》。

6.5.1 避免简单下溢

有些应用程序实际上执行了许多非常接近零的计算。这在计算残数或微分修正的算法中很常见。为获得在数值上安全的最大性能，需要采用扩展精度运算来执行关键计算。如果应用程序是单精度应用程序，可采用双精度执行关键计算。

示例：采用单精度的简单点积计算：

sum = 0

DO i = 1, n

sum = sum + a(i) * b(i)

END DO

如果 a(i) 和 b(i) 非常小，会出现很多下溢。通过强制计算采用双精度，计算点积时会具有更高的准确性，并且可避免出现下溢情况：

DOUBLE PRECISION sum