附录A IEEE浮点运算标准 在数值计算中,小数在内存中是以浮点数格式表示和参与运算的.浮点数是数字(或者说数 值)在内存中的一种存储格式,它和定点数是相对的。 A1浮点数与定点数 浮点数和定点数中的“点”指的是小数点 所谓定点数,就是指小数点的位置是固定的,不会向前或者向后移动.比如我们用4个字节 (32位字长)来存储无符号的定点数x(通常用4个字节存储单精度数),并且约定:前16位表示整 数部分,后16位表示小数部分,如下图所示: ·整数部分 ·小数部分 这种表示方法的优点是:整数部分和小数部分一目了然,非常直观 白对于整数,所有位都用来存储整数部分,所以一般采用定点数格式存储 但定点数格式有个很大的缺点,就是所能表示的数的取值范围比较小,比如在前面的例子中 (4个字节),所能表示的最大值和最小值(非)是 xm=1111111111111111.11111111111111112≈216=65536, xmin=0000000000000000.00000000000000012=2-16≈1.5×10-5 这在科学计算中显然是远远不够的.比如电子的质量大约是9×10-28克,用4位定点数格式就无 法表示.即使是采用8个字节(通常用8个字节存储双精度数),假定前32位表示整数,后32位表 示小数,则所能表示的的数的范围为(0除外) xm≈232≈4.3×109,Emn=2-32≈2.3×10-10. 为了克服这个缺点,人们发明了一种更加科学的存储格式,即浮点数格式,也就是通常所说的 科学计数法.该格式以指数形式存储数字,不但节省内存,也非常直观,而且所能表示的数的范围 也大大增加。 A.2EEE中的浮点数的表示方法 自计算机发明以来,曾出现许多中不同的浮点数表示方式,但目前最通用的是EEE二进制浮 点数算术标准(EEE Standard for Binary Floating-Point Arithmetic,.简称EEE754标准). IEEE7S4标准的主要起草者是加州大学伯克利分校数学系的William Kahan教授,他帮助Inte 公司设计了8087浮点处理器,并以此为基础形成了IEEE754标准,Kahan教授也因此获得了1987 年的图灵奖。 299
仅供课堂教学使用,请勿外传 附录 A IEEE 浮点运算标准 在数值计算中, 小数在内存中是以浮点数格式表示和参与运算的. 浮点数是数字(或者说数 值)在内存中的一种存储格式, 它和定点数是相对的. A.1 浮点数与定点数 浮点数和定点数中的“点”指的是小数点. 所谓定点数,就是指小数点的位置是固定的, 不会向前或者向后移动. 比如我们用 4 个字节 (32 位字长) 来存储无符号的定点数 x (通常用 4 个字节存储单精度数), 并且约定: 前 16 位表示整 数部分, 后 16 位表示小数部分, 如下图所示: 这种表示方法的优点是: 整数部分和小数部分一目了然, 非常直观. b 对于整数, 所有位都用来存储整数部分, 所以一般采用定点数格式存储. 但定点数格式有个很大的缺点, 就是所能表示的数的取值范围比较小. 比如在前面的例子中 (4 个字节), 所能表示的最大值和最小值 (非零) 是 xmax = 1111111111111111.11111111111111112 ≈ 2 16 = 65536, xmin = 0000000000000000.00000000000000012 = 2−16 ≈ 1.5 × 10−5 . 这在科学计算中显然是远远不够的. 比如电子的质量大约是 9 × 10−28 克, 用 4 位定点数格式就无 法表示. 即使是采用 8 个字节 (通常用 8 个字节存储双精度数), 假定前 32 位表示整数, 后 32 位表 示小数, 则所能表示的的数的范围为 (0 除外) xmax ≈ 2 32 ≈ 4.3 × 109 , xmin = 2−32 ≈ 2.3 × 10−10 . 为了克服这个缺点, 人们发明了一种更加科学的存储格式, 即浮点数格式, 也就是通常所说的 科学计数法. 该格式以指数形式存储数字, 不但节省内存, 也非常直观, 而且所能表示的数的范围 也大大增加. A.2 IEEE 中的浮点数的表示方法 自计算机发明以来, 曾出现许多中不同的浮点数表示方式, 但目前最通用的是 IEEE 二进制浮 点数算术标准 (IEEE Standard for Binary FloatingPoint Arithmetic, 简称 IEEE 754 标准). IEEE 754 标准的主要起草者是加州大学伯克利分校数学系的 William Kahan 教授, 他帮助 Intel 公司设计了 8087 浮点处理器, 并以此为基础形成了 IEEE 754 标准, Kahan 教授也因此获得了 1987 年的图灵奖. 299
,300 附录AEEE浮点运算标准 白William Kahan由于在浮点运算标准的制定上的杰出贡献,于1990年1月获得了图灵奖 通常一个浮点数由符号、尾数、基和指数组成,如: -0.3141592610×102,0.101012×23 这里要求小数点前面为零,小数点后面的数称为尾数.若尾数的首位数字不为0时,我们称其为正 规数(或规范化数),否则称为次正规数(或非规范化数).如0.31410×102是正规数而0.0031410× 10是次正规数.正规化表示方法可以使得每个浮点数的表示方式唯一,而且可以空出一个位置, 使得表示精度更高。 。EEE754标准中定义了表示浮点数的四种格式: ,两种基本的浮点数:单精度(32位字长)和双精度(64位字长). 其中单精度格式具有24位有效数字(二进制,而双精度格式具有53位有效数字(二进 制,相对于十进制来说,分别是7位(221≈10和16位(23≈1016)有效数字 两种扩展的浮点数:单精度扩展和双精度扩展。 EEE754标准中并未规定扩展格式的精度和大小,但它指定了最小精度和字长:单精度 扩展需43位字长以上,双精确度扩展需79位字长以上(64位有效数字).单精度扩展很 少使用,而对于双精确度扩展,不同的机器架构中有着不同的规定,有的为80位字长(如 X8O,有的为128位字长(如SPARC) 。一般来说,描述一个浮点数的三个基本要素为: ,基:计算机一般都以2为基; 尾数的位数:确定有效数字的位数,即精度 指数的位数:确定所能表示的数的范围 在EEE754标准中,浮点数是用二进制表示的,由三部分组成:符号(sgn,其值用s表示,指 数(exponent,其值用e表示)和尾数(fraction,,其值用f表示),见图AL.单精度数占32位字长 4个字节),第1位是符号位,第2至9位(8位字长)是指数位,最后23位是尾数.双精度数士 64位字长(8个字节),第1位是符号位,第2至12位(11位字长)是指数位,最后52位是尾数 IEEE单精度 1 8 23 智+ 数 (f) EEE双精度 1 11 52 图A1.EEE754中单精度格式与双精度格式的位模式 ·单精度格式:用8位字长的二进制数来表示指数,因此e的取值范围为0,25.当0≤e< 255时,按单精度格式存储的数,其对应的值是使用以下方法得到的: http://math.ecnu.edu.cn/-jypan
仅供课堂教学使用,请勿外传 · 300 · 附录 A IEEE 浮点运算标准 b William Kahan 由于在浮点运算标准的制定上的杰出贡献, 于 1990 年 1 月获得了图灵奖. 通常一个浮点数由符号、尾数、基和指数组成, 如: −0.3141592610 × 102 , 0.101012 × 2 3 . 这里要求小数点前面为零, 小数点后面的数称为尾数. 若尾数的首位数字不为 0 时, 我们称其为正 规数 (或规范化数), 否则称为次正规数 (或非规范化数). 如 0.31410 ×102 是正规数, 而 0.0031410 × 104 是次正规数. 正规化表示方法可以使得每个浮点数的表示方式唯一, 而且可以空出一个位置, 使得表示精度更高. • IEEE 754 标准中定义了表示浮点数的四种格式: 两种基本的浮点数: 单精度 (32 位字长) 和双精度 (64 位字长). 其中单精度格式具有 24 位有效数字 (二进制), 而双精度格式具有 53 位有效数字 (二进 制), 相对于十进制来说, 分别是 7 位 (2 24 ≈ 107 ) 和 16 位 (2 53 ≈ 1016) 有效数字. 两种扩展的浮点数: 单精度扩展和双精度扩展. IEEE 754 标准中并未规定扩展格式的精度和大小, 但它指定了最小精度和字长: 单精度 扩展需 43 位字长以上, 双精确度扩展需 79 位字长以上 (64 位有效数字). 单精度扩展很 少使用, 而对于双精确度扩展, 不同的机器架构中有着不同的规定, 有的为 80 位字长 (如 X86), 有的为 128 位字长 (如 SPARC). • 一般来说, 描述一个浮点数的三个基本要素为: 基: 计算机一般都以 2 为基; 尾数的位数: 确定有效数字的位数, 即精度; 指数的位数: 确定所能表示的数的范围. • 在 IEEE 754 标准中, 浮点数是用二进制表示的, 由三部分组成: 符号 (sign, 其值用 s 表示), 指 数 (exponent, 其值用 e 表示) 和尾数 (fraction, 其值用 f 表示), 见图 A.1. 单精度数占 32 位字长 (4 个字节), 第 1 位是符号位, 第 2 至 9 位 (8 位字长) 是指数位, 最后 23 位是尾数. 双精度数占 64 位字长 (8 个字节), 第 1 位是符号位, 第 2 至 12 位 (11 位字长) 是指数位, 最后 52 位是尾数. 图 A.1. IEEE 754 中单精度格式与双精度格式的 位模式 • 单精度格式: 用 8 位字长的二进制数来表示指数, 因此 e 的取值范围为 [0, 255]. 当 0 ≤ e < 255 时, 按单精度格式存储的数, 其对应的值是使用以下方法得到的: http://math.ecnu.edu.cn/~jypan
A.2IEEE中的浮点数的表示方法 ·301 将二进制基数点(小数点)插入到尾数∫最高有效位的左侧,并将一个隐含位插人到二进 制基数点的左侧,从而得到的是一个二进制带分数(整数加小数) 由此构成的带分数就是单精度格式有效数字.隐含位的值并不是显式指定的(即不存储,而 是通过指数e的值来隐式指定: ·当0<e<255时.表示该数为一进制正规数.此时隐含位设为1. ·当e=0时,表示该数为二进制次正规数,隐含位设为0. 凸由于引入了隐含位(伪了尽可能地增加所能表示的数的精度),这里的正规数概念与前 面的定义有点区别,因此我们加上“二进制”三个字 单精度格式位模式中的尾数只有23位,但由于使用了隐含位,所以能提供24位有效数字(仁 进制).因此,在EEE中,单精度数的表示方法为 (-1)°×1.f×2-127(仁进制正规数) (-1)°×0.f×2-126仁进制次正规数 完整的对应关系是 单精度格式位模式 值 0<e<255 (-1)°×1.f×2e-127(二进制正规数物 e=0,f≠0 (-1)×0.f×2-126(仁进制次正规数 e=0,f=0 (-1)×0.0(有符号的零) e=255,f=0,8=0 +inf(正无穷大) e=255,f=0,8=1 inf(负无穷大) e=255,f≠0 NaN(非数、非确定值 其中127是单精度格式的指数偏移值(exponent bias),在IEEE标准中,这个值定义为2(指数位长-1) 1.所以对于单精度格式,指数偏移值就是28-1-1=127,而对于双精度格式,这个值为 21- -1=1023. ·双精度格式:与单精度格式类似,完整的对应关系是 双精度格式位模式 0<e<2047 (-1)×1.f×2-102(仁进制正规数) e=0,f≠0 (-1)°×0.f×21022(仁进制次正规数) e=0,f=0 (-1)°×0.0(有符号的零) e=2047,f=0,8=0 +inf(正无穷大) e=2047,f=0,s=1 inf(负无穷大 e=2047,f≠0 NaN(非数、非确定值) http://math.ecnu.edu.cn/-jypan
仅供课堂教学使用,请勿外传 A.2 IEEE 中的浮点数的表示方法 · 301 · 将二进制基数点 (小数点) 插入到尾数 f 最高有效位的左侧, 并将一个隐含位插入到二进 制基数点的左侧, 从而得到的是一个二进制带分数 (整数加小数). 由此构成的带分数就是单精度格式有效数字. 隐含位的值并不是显式指定的 (即不存储), 而 是通过指数 e 的值来隐式指定: • 当 0 < e < 255 时, 表示该数为二进制正规数, 此时隐含位设为 1. • 当 e = 0 时, 表示该数为二进制次正规数, 隐含位设为 0. b 由于引入了隐含位 (为了尽可能地增加所能表示的数的精度), 这里的正规数概念与前 面的定义有点区别, 因此我们加上 “二进制” 三个字. 单精度格式位模式中的尾数只有 23 位, 但由于使用了隐含位, 所以能提供 24 位有效数字 (二 进制). 因此, 在 IEEE 中, 单精度数的表示方法为 (−1)s × 1.f × 2 e−127 (二进制正规数) (−1)s × 0.f × 2 −126 (二进制次正规数) 完整的对应关系是 单精度格式位模式 值 0 < e < 255 (−1)s × 1.f × 2 e−127 (二进制正规数) e = 0, f ̸= 0 (−1)s × 0.f × 2 −126 (二进制次正规数) e = 0, f = 0 (−1)s × 0.0 (有符号的零) e = 255, f = 0, s = 0 +inf (正无穷大) e = 255, f = 0, s = 1 inf (负无穷大) e = 255, f ̸= 0 NaN (非数、非确定值) 其中 127 是单精度格式的指数偏移值 (exponent bias), 在IEEE标准中, 这个值定义为2 (指数位长−1)− 1. 所以对于单精度格式, 指数偏移值就是 2 8−1 − 1 = 127, 而对于双精度格式, 这个值为 2 11−1 − 1 = 1023. • 双精度格式: 与单精度格式类似, 完整的对应关系是 双精度格式位模式 值 0 < e < 2047 (−1)s × 1.f × 2 e−1023 (二进制正规数) e = 0, f ̸= 0 (−1)s × 0.f × 2 −1022 (二进制次正规数) e = 0, f = 0 (−1)s × 0.0 (有符号的零) e = 2047, f = 0, s = 0 +inf (正无穷大) e = 2047, f = 0, s = 1 inf (负无穷大) e = 2047, f ̸= 0 NaN (非数、非确定值) http://math.ecnu.edu.cn/~jypan
·302: 附录A IEEE浮点运算标准 例A1单精度格式所能表示的十进制数范围 ·最大二进制正规数为7F7 FFFFF16=3.40282347×1038 ·最小正的二进制正规数为0800016=117549435×10- ·最大二进制次正规数为007 FFFFF16=1.17549421×10- ·最小(正的)二进制次正规数为0000000116=1.40129846×10-45 由此可见,浮点数所能表示的数的范围比定点数要大很多, 例A2双精度格式所能表示的十进制数范围。 ·最大二进制正规数为7 FEFFFFF FFFFFFFF16=1.7976931348623157×100 。最小(正的)二进制正规数为001000000000000016=2.2250738585072014×10-308 ·最大二进制次正规数为000 FFFFF FFFFFFFF16=2.2250738585072009×10-08 ·最小(正的)二进制次正规数为000000000000000116=4.9406564584124654×10-2 ·3FF000000000000016=1(补码) 白在MATLAB中,hex2num可以将一个由16个16进制数组成的字符串转化为其所对应 的浮点数(根据IEEE标准,类似的命令有hex2dec,bin2dec,base2dec,num2hex, dec2bin,... 例A.3把二进制数(1001.0101)2转换成十进制数. (1001.0101)2=1×23+0×22+0×21+1×20 +0×2-1+1×2-2+0×2-3+1×2-2 =9.312510 例A4把十进制数13.12510转换成二进制数 整数部分:1310=11012(银转相除法) 小数部分: 。0.125×2=0.25,整数位是0→.0: ·0.25×2=0.5,整数位是0+.00 ·0.5×2=1,整数位是1→.001;(小数部分已变为0,运算结束) 所以13.12510=1101.0012 一个十进制数能否用二进制浮点数精确表示,关键在于小数部分】 例A.5十进制数0.110能否用二进制数精确表示? 。0.1×2=0.2,整数位是0→.0: ·0.2×2=0.4,整数位是0→.00: http://math.ecnu.edu.cn/-jypan
仅供课堂教学使用,请勿外传 · 302 · 附录 A IEEE 浮点运算标准 例 A.1 单精度格式所能表示的十进制数范围. • 最大二进制正规数为 7F7FFFFF16 = 3.40282347 × 1038 • 最小 (正的) 二进制正规数为 0080000016 = 1.17549435 × 10−38 • 最大二进制次正规数为 007FFFFF16 = 1.17549421 × 10−38 • 最小 (正的) 二进制次正规数为 0000000116 = 1.40129846 × 10−45 由此可见, 浮点数所能表示的数的范围比定点数要大很多. 例 A.2 双精度格式所能表示的十进制数范围. • 最大二进制正规数为 7FEFFFFF FFFFFFFF16 = 1.7976931348623157 × 10308 • 最小 (正的) 二进制正规数为 00100000 0000000016 = 2.2250738585072014 × 10−308 • 最大二进制次正规数为 000FFFFF FFFFFFFF16 = 2.2250738585072009 × 10−308 • 最小 (正的) 二进制次正规数为 00000000 0000000116 = 4.9406564584124654 × 10−324 • 3FF00000 0000000016 = 1 (补码) b 在 MATLAB 中, hex2num 可以将一个由 16 个 16 进制数组成的字符串转化为其所对应 的浮点数 (根据 IEEE 标准), 类似的命令有 hex2dec, bin2dec, base2dec, num2hex, dec2bin, ... 例 A.3 把二进制数 (1001.0101)2 转换成十进制数. (1001.0101)2 =1 × 2 3 + 0 × 2 2 + 0 × 2 1 + 1 × 2 0 + 0 × 2 −1 + 1 × 2 −2 + 0 × 2 −3 + 1 × 2 −2 =9.312510 例 A.4 把十进制数 13.12510 转换成二进制数. 整数部分: 1310 = 11012 (辗转相除法) 小数部分: • 0.125 × 2 = 0.25, 整数位是 0 → .0; • 0.25 × 2 = 0.5, 整数位是 0 → .00; • 0.5 × 2 = 1, 整数位是 1 → .001; (小数部分已变为 0, 运算结束) 所以 13.12510 = 1101.0012 一个十进制数能否用二进制浮点数精确表示, 关键在于小数部分. 例 A.5 十进制数 0.110 能否用二进制数精确表示? • 0.1 × 2 = 0.2, 整数位是 0 → .0; • 0.2 × 2 = 0.4, 整数位是 0 → .00; http://math.ecnu.edu.cn/~jypan
A3IEEE中的浮点数运算 303 ·0.4×2=0.8,整数位是0→.000 ·0.8×2=1.6,整数位是1→.0001: ·0.6×2=1.2,整数位是1→.00011 ·0.2×2=0.4,整数位是0→.000110: 得到一个无限循环的二进制小数,显然用有限位字长是无法表示的,因此0.10无法用EEE754 浮点数精确表示. 同理可知0.2.0.4.0.6.0.8.0.3.0.7.0.9也是无法精确表示的.故0.1至0.9的9个小数中.只 有0.5可以精确表示. 例A6能用二进制数精确表示的十进制小数.易知 0.12=2=0.5 0.012=28=0.25 0.0012=20=0.125 0.00012=20=0.0625 0.000012=20=0.03125 0.0000012=20=0.015625 0.00000012=20=0.0078125 0.000000012=20=0.00390625 由此可知,一个十进制小数要能用浮点数精确表示,最后一位必须是5.当然这是必要条件,并非 充分条件.如0.35就无法精确表示. 例A7N位二进制小数能精确表示的非零十进制小数总共有多少个? ·1位二进制小数能精确表示的有20=1个(0.12=0.510: ·2位二进制小数能精确表示的有21=2个(0.012=0.2510,0.112=0.7510: ·3位二进制小数能精确表示的有2=4个 ·N位二进制小数能精确表示的有2V-1个 所以N位二进制小数能精确表示的十进制小数总共有2V-1个 A.3IEEE中的浮点数运算 。EEE754标准也定义了浮点数的运算规侧: -加.减、乘.除。平方根.余数、将浮点格式的数含人为整数值、在不同浮点格式之向 转换、在浮点和整数格式之间转换以及比较:EEE对以上浮点运算的准确度作了规定: http://math.ecnu.edu.cn/-jypan
仅供课堂教学使用,请勿外传 A.3 IEEE 中的浮点数运算 · 303 · • 0.4 × 2 = 0.8, 整数位是 0 → .000; • 0.8 × 2 = 1.6, 整数位是 1 → .0001; • 0.6 × 2 = 1.2, 整数位是 1 → .00011; • 0.2 × 2 = 0.4, 整数位是 0 → .000110; • · · · · · · 得到一个无限循环的二进制小数, 显然用有限位字长是无法表示的, 因此 0.110 无法用 IEEE 754 浮点数精确表示. 同理可知 0.2, 0.4, 0.6, 0.8, 0.3, 0.7, 0.9 也是无法精确表示的. 故 0.1 至 0.9 的 9 个小数中, 只 有 0.5 可以精确表示. 例 A.6 能用二进制数精确表示的十进制小数. 易知 0.12 = 2−1 10 = 0.5 0.012 = 2−2 10 = 0.25 0.0012 = 2−3 10 = 0.125 0.00012 = 2−4 10 = 0.0625 0.000012 = 2−5 10 = 0.03125 0.0000012 = 2−6 10 = 0.015625 0.00000012 = 2−7 10 = 0.0078125 0.000000012 = 2−8 10 = 0.00390625 · · · · · · 由此可知, 一个十进制小数要能用浮点数精确表示, 最后一位必须是 5. 当然这是必要条件, 并非 充分条件. 如 0.35 就无法精确表示. 例 A.7 N 位二进制小数能精确表示的非零十进制小数总共有多少个? • 1 位二进制小数能精确表示的有 2 0 = 1 个 (0.12 = 0.510); • 2 位二进制小数能精确表示的有 2 1 = 2 个 (0.012 = 0.2510, 0.112 = 0.7510); • 3 位二进制小数能精确表示的有 2 2 = 4 个 • · · · · · · • N 位二进制小数能精确表示的有 2 N−1 个 所以 N 位二进制小数能精确表示的十进制小数总共有 2 N−1 个. A.3 IEEE 中的浮点数运算 • IEEE 754 标准也定义了浮点数的运算规则: 加、减、乘、除、平方根、余数、将浮点格式的数舍入为整数值、在不同浮点格式之间 转换、在浮点和整数格式之间转换以及比较: IEEE 对以上浮点运算的准确度作了规定: http://math.ecnu.edu.cn/~jypan
,304 附录AEEE浮点运算标准 求余和比较运算必须精确无误.其他运算必须向其目标提供精确的结果,除非没有此类 结果,或者该结果不满足目标格式,此时运算必须按照下面介绍的舍入棋式对精确结果 进行最低限度的修改,并将经过修改的结果提供给运算的目标。 ·在十进制字符串和两种基本浮点格式之一的二进制浮点数之间进行转换的准确度、单 性和一致性要求:对于在指定范围内的操作数,这些转换必须生成精确的结果(如果可 能的话),或者按照规定的舍入模式,对此类精确结果进行最低限度的修改.对于不在指 定范围内的操作数,这些转换生成的结果与精确结果之间的差值不得超过取决于舍人模 式的指定误差. 五种类型的EEE浮点异常:无效运算(如0/0,/∞等),被零除,上溢,下溢和不精确 以及用于向用户指示发生这些类型异常的条件 四种舍入模式:设x是所要表示的数, ()就近舍人:用最接近x的可表示的值来代替,类似于整数的四舍五入.如果x正 好在两个相邻的可表示值的中间,则首选二进制“偶数”(二进制最后一位为0: (2)向下舍人:用不大于x的可表示的值来代替(向负无穷大方向截断): (3)向上舍人:用不小于x的可表示的值来代替(向正无穷大方向截断) (④)向0舍人:当x>0时采用向下舍入,当x<0时采用向上舍入 凸我们将后面三种舍入模式统称为截断。 凸不同编译器对舍入可能有不同的处理方式, 。下溢 当运算结果非常小时,就会发生下溢.下表是下溢值 目标的精度 下溢阀值 单精度 最小正规数1.17549435×10-3 最大次正规数1.17549421×10-3然 双精度 最小正规数2.250738585072014×10- 最大次正规数2.2250738585072009×10-08 IEEE算法处理下溢的方式是渐进下溢:当生成的正确结果的数量级低于最小正正规数时,就 会生成次正规数,而不是返回零 ·机器精度:将l.0与大于1.0的最小浮点数之间的距离记为em,它的一半称为unit roundo乐, 记为,它是计算机表示一个浮点数时的相对误差界, )=1+或四=1+间≤ 这里A(x)表示x在计算机中实际存储的IEEE浮点数 在EEE标准下,单精度和双精度浮点运算的最大相对误差©u分别为 http://math.ecnu.edu.cn/-jypan
仅供课堂教学使用,请勿外传 · 304 · 附录 A IEEE 浮点运算标准 求余和比较运算必须精确无误. 其他运算必须向其目标提供精确的结果, 除非没有此类 结果, 或者该结果不满足目标格式, 此时运算必须按照下面介绍的舍入模式对精确结果 进行最低限度的修改, 并将经过修改的结果提供给运算的目标. 在十进制字符串和两种基本浮点格式之一的二进制浮点数之间进行转换的准确度、单 一性和一致性要求: 对于在指定范围内的操作数, 这些转换必须生成精确的结果 (如果可 能的话), 或者按照规定的舍入模式, 对此类精确结果进行最低限度的修改. 对于不在指 定范围内的操作数, 这些转换生成的结果与精确结果之间的差值不得超过取决于舍入模 式的指定误差. 五种类型的 IEEE 浮点异常: 无效运算 (如 0/0, ∞/∞ 等), 被零除, 上溢, 下溢和不精确, 以及用于向用户指示发生这些类型异常的条件. 四种舍入模式: 设 x 是所要表示的数, (1) 就近舍入 : 用最接近 x 的可表示的值来代替, 类似于整数的四舍五入. 如果 x 正 好在两个相邻的可表示值的中间, 则首选二进制 “偶数” (二进制最后一位为 0); (2) 向下舍入 : 用不大于 x 的可表示的值来代替 (向负无穷大方向截断); (3) 向上舍入 : 用不小于 x 的可表示的值来代替 (向正无穷大方向截断); (4) 向 0 舍入 : 当 x > 0 时采用向下舍入, 当 x < 0 时采用向上舍入. b 我们将后面三种舍入模式统称为 截断. b 不同编译器对舍入可能有不同的处理方式. • 下溢 当运算结果非常小时, 就会发生下溢. 下表是下溢阈值. 目标的精度 下溢阈值 单精度 最小正规数 1.17549435 × 10−38 最大次正规数 1.17549421 × 10−38 双精度 最小正规数 2.2250738585072014 × 10−308 最大次正规数 2.2250738585072009 × 10−308 IEEE 算法处理下溢的方式是渐进下溢: 当生成的正确结果的数量级低于最小正正规数时, 就 会生成次正规数, 而不是返回零. • 机器精度: 将 1.0 与大于 1.0 的最小浮点数之间的距离记为 εm. 它的一半称为 unit roundoff, 记为 εu, 它是计算机表示一个浮点数时的相对误差界, fl(x) = x(1 + δ) 或 fl(x) = x 1 + δ , |δ| ≤ εu. 这里 fl(x) 表示 x 在计算机中实际存储的 IEEE 浮点数. 在 IEEE 标准下, 单精度和双精度浮点运算的最大相对误差 εu 分别为 http://math.ecnu.edu.cn/~jypan
A.4浮点运算舍人误差分析 ·305 精度 最大相对误差 单精度 2-24≈5.960464×10-8 双精度 2-53≈1.110223×10-16 凸如果采用的不是就近舍入模式,而是其他三种合入模式(即截断),则最大相对误差为 Em 么有的文献中称cu为机器精度(machine epsilon,machine precision,or macheps),如 Demmel[30,LAPACK.Scilab,Wikipedia.也有的文献称em为机器精度,如Higham[68, MATLAB,.Mathematica.我们采用前面一种方式,即“机器精度”指的是eau. 例A.8假定要使用只有三个精度位的二进制算法.那么,最大相对误差为23.在任意两个 2的幂之间,只有23-1=7个可表示数字,如下图所示 01212 16 2.12021 02 数轴显示了数字之间的差距是随着指数增加而加倍增加的 在正EE单精度格式中,两个最小正次正规数之间的差大约是105,而两个最大有限数之间 的数量级差大约是1031! 凸精确是偶然的,误差是必然的.做数值算法,惟一能做的就是尽量使误差的传播和累积能够 得到有效的控制 A4浮点运算合人误差分析 由于计算机无法精确表示所有的浮点数,在做浮点运算时,如果计算结果无法精确表示,此时 就会产生的误差,这就是浮点运算的含人误差.根据EEE浮点运算标准,如果a⊙b的结果无法精 确表示,则用一个最接近的浮点数来代替(浮点运算时一般采用就近舍入模式),记为(a⊙).这 里的⊙表示加、减、乘、除四种运算符. 在不考虑溢出的情况下,我们有 laob=aob1+)或ao)=i aob 其中6表示浮点运算的相对误差,满足1≤©,公式(A1)是分析浮点运算舍人误差的基础(标 准模型)[68,Pg401.EEE浮点运算标准同时也指出,对于开根号运算,产生的误差同样也满足 http://math.ecnu.edu.cn/-jypan
仅供课堂教学使用,请勿外传 A.4 浮点运算舍入误差分析 · 305 · 精度 最大相对误差 单精度 2 −24 ≈ 5.960464 × 10−8 双精度 2 −53 ≈ 1.110223 × 10−16 b 如果采用的不是就近舍入模式, 而是其他三种舍入模式 (即截断), 则最大相对误差为 εm. b 有的文献中称 εu 为机器精度 (machine epsilon , machine precision , or macheps ), 如 Demmel [30], LAPACK, Scilab, Wikipedia. 也有的文献称 εm 为机器精度, 如 Higham [68], MATLAB, Mathematica. 我们采用前面一种方式, 即 “机器精度” 指的是 εu. 例 A.8 假定要使用只有三个精度位的二进制算法. 那么, 最大相对误差为 2 −3 . 在任意两个 2 的幂之间, 只有 2 3 − 1 = 7 个可表示数字, 如下图所示. 数轴显示了数字之间的差距是随着指数增加而加倍增加的. 在 IEEE 单精度格式中, 两个最小正次正规数之间的差大约是 10−45 , 而两个最大有限数之间 的数量级差大约是 1031! b 精确是偶然的, 误差是必然的. 做数值算法, 惟一能做的就是尽量使误差的传播和累积能够 得到有效的控制. A.4 浮点运算舍入误差分析 由于计算机无法精确表示所有的浮点数, 在做浮点运算时, 如果计算结果无法精确表示, 此时 就会产生的误差, 这就是浮点运算的舍入误差. 根据 IEEE 浮点运算标准, 如果 a ⊙ b 的结果无法精 确表示, 则用一个最接近的浮点数来代替 (浮点运算时一般采用就近舍入模式), 记为 fl(a ⊙ b). 这 里的 ⊙ 表示加、减、乘、除四种运算符. 在不考虑溢出的情况下, 我们有 fl(a ⊙ b) = (a ⊙ b)(1 + δ) 或 fl(a ⊙ b) = a ⊙ b 1 + δ , (A.1) 其中 δ 表示浮点运算的相对误差, 满足 |δ| ≤ εu. 公式 (A.1) 是分析浮点运算舍入误差的基础 (标 准模型) [68, page 40]. IEEE 浮点运算标准同时也指出, 对于开根号运算, 产生的误差同样也满足 http://math.ecnu.edu.cn/~jypan
·306. 附录A IEEE浮点运算标准 一个比较有趣的问题是何时两个浮点数能进行精确的相减运算。 定理A1(Ferguson,,1995)[68设浮点数x和y满足 e(z-y)s minfe(z),e(y)} (A2 其中()表示一个正规浮点数的指数则(红-)=工一.(假定一y不会下溢) 白这里只考虑浮点运算误差,所以假定x,y都是计算机可以精确表示的.事实上,由(A2)可 知,x和y的指数至多差1. 上面的定理对任何基都成立.如果基是2,则有下面的结论 推论A.2(Sterbenz,197④设浮,点数x和y满足 /2≤x≤2y, 则(红一)=x一弘.(假定x一y不会下溢) 下面的结论对估计浮点运算的舍人误差非常有用[68,pg©63), 引理A3设耐≤eu且neu<1,则 ia+sr=1+,a≤a色 其中:=士1 例A)(多项式运算的金人误差分析:已知多项式)=三对于给定的五分析计算 p(c)的值时的舍入误差, 解.在对多项式p叫)=anx+an-1x-1…+a1x+o求值时,我们通常采用Horner法则 算法Al.Horner法则 1:p=a 2:fori=n-1:-1:0 do p=x+p+ai 4:end for 先看一个具体的例子.设p(x)=(任-2),观测工=2附近的舍入误差.我们通过画图来显示 误差情祝(见图A.2.观察发现,用Hoer法则求值时会在x=2附近会出现“噪声带 下面我们基于舍入误差分析模型(A.1)来分析多项式求值的误差情况.带舍人误差的Hore 算法可写为 http://math.ecnu.edu.cn/-jypan
仅供课堂教学使用,请勿外传 · 306 · 附录 A IEEE 浮点运算标准 (A.1). 一个比较有趣的问题是何时两个浮点数能进行精确的相减运算. 定理 A.1 (Ferguson, 1995) [68] 设浮点数 x 和 y 满足 e(x − y) ≤ min{e(x), e(y)} (A.2) 其中 e(·) 表示一个正规浮点数的指数. 则 fl(x − y) = x − y. (假定 x − y 不会下溢) b 这里只考虑浮点运算误差, 所以假定 x, y 都是计算机可以精确表示的. 事实上, 由 (A.2) 可 知, x 和 y 的指数至多差 1. 上面的定理对任何基都成立. 如果基是 2, 则有下面的结论. 推论 A.2 (Sterbenz, 1974) 设浮点数 x 和 y 满足 y/2 ≤ x ≤ 2y, 则 fl(x − y) = x − y. (假定 x − y 不会下溢) 下面的结论对估计浮点运算的舍入误差非常有用 [68, page 63]. 引理 A.3 设 |δi | ≤ εu 且 nεu < 1, 则 Yn i=1 (1 + δi) ρi = 1 + θn, |θn| ≤ γn ≜ nεu 1 − nεu , 其中 ρi = ±1. 例 A.9 (多项式运算的舍入误差分析) : 已知多项式 p(x) = Pn k=0 akx k . 对于给定的 x, 分析计算 p(x) 的值时的舍入误差. 解. 在对多项式 p(x) = anx n + an−1x n−1 · · · + a1x + a0 求值时, 我们通常采用 Horner 法则. 算法 A.1. Horner 法则 1: p = an 2: for i = n − 1 : −1 : 0 do 3: p = x ∗ p + ai 4: end for 先看一个具体的例子. 设 p(x) = (x − 2)9 , 观测 x = 2 附近的舍入误差. 我们通过画图来显示 误差情况 (见图 A.2). 观察发现, 用 Horner 法则求值时会在 x = 2 附近会出现 “噪声带”. 下面我们基于舍入误差分析模型 (A.1) 来分析多项式求值的误差情况. 带舍入误差的 Horner 算法可写为 http://math.ecnu.edu.cn/~jypan
A4浮点运算舍人误差分析 一 图A2.分别利用Horner 法则(左图)和y=(-2)9(右图)在区间1.92,2.08)上的 8000个等分点上求值的结果 算法A2.带入误差的Horner法则 1:p=On 2:fori=n-1:-1:0 do p=(*p)(1+)+a(1+)%≤e,1l≤e 4:end for 所以得到的最终计算结果为 p()= 假定j·ug1,则有 (1+d1).(1+6)≤(1+em≤1+j·eu (1+61)…(1+6)≥(1-em≥1-j·eu 于是(A.3)可改写为 a=∑(1+)ara,其中ls2n A到 所以()可看作是另一个多项式的精确值,且这个多项式的系数是原多项式系数的一个小的 扰动.这说明多项式计算的Horner算法是向后稳定的,且系数的向后相对误差不超过2n·ew,因 此,绝对误差为 Ip(x》-p(xl= a-apr n. http://math.ecnu.edu.cn/-jypan
仅供课堂教学使用,请勿外传 A.4 浮点运算舍入误差分析 · 307 · 图 A.2. 分别利用 Horner 法则(左图)和 y = (x − 2)9(右图)在区间 [1.92, 2.08] 上的 8000 个等分点上求值的结果 算法 A.2. 带舍入误差的 Horner 法则 1: p = an 2: for i = n − 1 : −1 : 0 do 3: p = ((x ∗ p)(1 + δi) + ai)(1 + ˜δi) % |δi | ≤ εu, | ˜δi | ≤ εu 4: end for 所以得到的最终计算结果为 fl(p(x)) = nX−1 i=0 (1 + ˜δi) Y i−1 j=0 (1 + δj )(1 + ˜δj ) aix i + nY−1 j=0 (1 + δj )(1 + ˜δj ) anx n . (A.3) 假定 j · εu ≪ 1, 则有 (1 + δ1)· · ·(1 + δj ) ≤ (1 + εu) j ≲ 1 + j · εu, (1 + δ1)· · ·(1 + δj ) ≥ (1 − εu) j ≥ 1 − j · εu. 于是 (A.3) 可改写为 fl(p(x)) = Xn i=0 1 + ˆδi aix i ≜ Xn i=0 a˜ix i , 其中 | ˆδi | ≤ 2n · εu. (A.4) 所以 fl(p(x)) 可看作是另一个多项式的精确值, 且这个多项式的系数是原多项式系数的一个小的 扰动. 这说明多项式计算的 Horner 算法是向后稳定的, 且系数的向后相对误差不超过 2n · εu. 因 此, 绝对误差为 | fl(p(x)) − p(x)| = Xn i=0 a˜ix i − Xn i=0 aix i = Xn i=0 ˆδiaix i ≤ Xn i=0 ˆδiaix i ≤ 2n · εu Xn i=0 |aix i |. http://math.ecnu.edu.cn/~jypan
·308. 附录A IEEE浮点运算标准 相对误差为 lA(p(z))-p(z)l Ip() ≤2m A5课后习题 蛛习11考虑抽样数据1,2,其均值为元=1乃,样本方差为无偏估计 试证明 2-4-m 对于数值计算,上述两种计算方法哪一个更可靠,为什么: 练习12根据浮点运算的舍入误差模型(A1),在不考虑溢出的情况下,证明[68) http://math.ecnu.edu.cn/-jypan
仅供课堂教学使用,请勿外传 · 308 · 附录 A IEEE 浮点运算标准 相对误差为 | fl(p(x)) − p(x)| |p(x)| ≤ 2n · εu Pn i=0 |aix i | Pn i=0 aix i . □ A.5 课后习题 练习 1.1 考虑抽样数据 x1, x2, . . . , xn, 其均值为 x¯ = 1 n Xn i=1 xi , 样本方差为 (无偏估计) σ 2 = 1 n − 1 Xn i=1 (xi − x¯) 2 . 试证明: σ 2 = 1 n − 1 Xn i=1 (x 2 i − nx¯ 2 ). 对于数值计算, 上述两种计算方法哪一个更可靠, 为什么? 练习 1.2 根据浮点运算的舍入误差模型 (A.1), 在不考虑溢出的情况下, 证明 [68] fl Xn i=1 xiyi ! = Xn i=1 xiyi(1 + θi), 其中 |θi | ≤ γi ≜ iεu 1 − iεu , i = 1, 2, . . . , n. http://math.ecnu.edu.cn/~jypan