无论是解决素数分布问题的基础科学、实现全球天气预 的应用科学、验证导弹飞行实验的工程项目,还是日常生活中一张图片的显示(依赖于余弦变换)、一个通信的建立(依赖于语音采样)都离不开浮点计算,而对于计算而言,最基本的要求就是其结果的绝对正确性,但这个基本要求在计算机系统中很难得到保证.目前,科学计算大多使用浮点硬件,且都遵循IEEE-754浮点标准[,此标准定义的浮点数在数轴上呈离散态,在计算过程中需要对计算结果进行近似表示,而计算总是伴随着舍入,增加了浮点不精确甚至异常的可能.例如:由于浮点数转整数出现的整数溢出异常,欧洲Ariane 5火箭在1996年发射时出现了严重的升空自爆现象[,造成了巨额的经济损失.2010年,由于控制软件异常导致的丰田汽车刹车问题,同样给个人带来了极大的安全隐患[.这些事件都表明了浮点数计算异常的危害性及其难以控制性.此外,即使没有浮点计算异常,在常规计算中也可能出现非常规情况,即所谓的病态条件[,如文献[中的函数:
$f(x) = frac{{tan (sin (x)) – sin (tan (x))}}{{{x^7}}}.$
该函数在x=0.0处的真值为1/30,而利用4种不同计算软件计算的结果却存在较大差异,如图1所示.
图1(Fig.1)
图2
浮点数溢出
Fig.2
Floating-Point overflow
需要注意的是,下溢不是一个严重问题,通常看作为机器0.
DZE异常为当除数为0而被除数为有限数时触发的异常,也泛指有限数运算导致无穷结果的异常,例如4.0/0.0,log(0,0)等.在浮点函数实际运算过程中,一般只被除法指令FDIVS/FDIVD触发.
INE异常为由操作数或计算结果的不精确表示触发的异常,不精确表示的主要原因是浮点格式的有效位数不足以容纳精确的计算结果,如除不尽、两个浮点数相乘等情况都会触发该异常.由于INE异常会在浮点运算中经常被触发,在一般情况下总是被屏蔽.
若当前浮点操作的操作数中存在非规格化数,则触发DNO异常.
关于DNO异常的处理可以限定为以下情况:
1) 具有非规格化操作数的浮点指令不会产生上溢、下溢或非精确结果自陷;
2) 除以一个非规格化数根据不同情况可以作为DZE异常或INV异常;
3) 非规格化数乘以无穷大作为INV异常;
4) 对一个负的非规格化数求平方根将产生一个负0结果,而不产生异常;
5) 非规格化数被当作0,不会引起DNO的自陷;
6) 在浮点函数实际运算过程中,非规格化数多数被当作0.
上述6类异常都将会被最基本的浮点操作(加、减、乘、除)或浮点函数(sqrt,sin,pow,fp_class等)所触发.首先定义一个函数E:经过某一个运算(基本浮点操作或者浮点函数)后,可能出现的状态及其触发条件.对于浮点数OP1和OP2,经过E后有如下的规则:
$begin{array}{l}
E(OP1Theta OP2) = left{ {begin{array}{*{20}{l}}
{Invalid,{rm{ if }}(OP1{rm{ or }}OP2) = = NaN{rm{ or }}inf }\
{Overflow,{rm{ if }}|OP1Theta OP2| > {rm{MAX}}}\
{Underflow,{rm{ if }}0 < |OP1Theta OP2| < {rm{MIN}}}\
{Denormal,{rm{ if }}0 < |OP1| < {rm{MIN or }}0 < |OP2| < {rm{MIN}}}\
{OP1Theta OP2,{rm{ otherwise}}}
end{array}} right.,Theta = { + , – , times } ,\
E(OP1/OP2) = left{ {begin{array}{*{20}{l}}
{Invalid,{rm{ if }}((OP1{rm{ or }}OP2) = = NaN{rm{ or }}inf ){rm{ or }}(OP1 = 0{rm{ and }}OP2 = 0)}\
{Overflow,{rm{ if }}|OP1| > left| {OP2} right| cdot {rm{MAX}}}\
{Underflow,{rm{ if }}0 < |OP1| < |OP2| cdot {rm{MIN}}}\
{Denormal,{rm{ if }}0 < |OP1| < {rm{MIN or }}0 < |OP2| < {rm{MIN}}}\
{Divide{rm{ }}by{rm{ }}Zero,{rm{ if }}(OP1 ne 0{rm{ or }}NaN{rm{ or }}inf ){rm{ and }}OP2 = 0}\
{OP1/OP2,{rm{ otherwise}}}
end{array}} right..
end{array}$
四则运算是浮点计算的基础运算,了解经过四则运算可能产生的各种异常,是后续对浮点函数异常分析的基础.虽然每个函数由于自身的功能不同,可能导致的异常有所不同,但回到本源都是由基本运算引起的异常.部分浮点函数异常状况如下:
$begin{array}{l}
E(sqrt(OP1)) = left{ {begin{array}{*{20}{l}}
{Invalid,{rm{ if }}(OP1 < 0){rm{ or }}(OP1 = = NaN{rm{ or }}inf )}\
{z,{rm{ and }}z cdot z = x,{rm{ otherwise}}}
end{array}} right.,\
E(sin (OP1)) = left{ {begin{array}{*{20}{l}}
{Invalid,{rm{ if }}OP1 = = NaN{rm{ or }}inf }\
{z,{rm{ otherwise}}}
end{array}} right.,\
E(cot (OP1)) = left{ {begin{array}{*{20}{l}}
{Invalid,{rm{ if }}OP1 = = NaN{rm{ or }}inf }\
{Divide{rm{ }}by{rm{ }}Zero,{rm{ if }}OP1 = 0}\
{Overflow,{rm{ if }}|OP1| le 0x0003000000000000}\
{z,{rm{ otherwise}}}
end{array}} right.,\
E(fp_class(OP1)) = z.
end{array}$
对于上述规则,函数和浮点操作可以相互嵌套,如E(sin(OP1/OP2)),E(sqrt(OP1)Qcot(OP1))等.无论是以何种方式嵌套出现,整体的E都是各子调用的并集,即:
E(sqrt(OP1)Qcot(OP1))=E(sqrt(OP1))èE(cot(OP1))èE(OP1QOP2).
2. 分段式异常处理
该处理流程在浮点函数实现之前完成对浮点异常的编码,将每一类异常编码为可表示的浮点数;在此基础上对输入参数进行异常检测,即,处理可能出现的INV异常;在正常运算的过程中,保持对浮点操作特殊指令的检测,将处理重点集中在INF异常和DZE异常;最后,对计算结果进行检测,处理可能存在的FPF异常和DNO异常.其具体算法见算法1.
图3(Fig.3)
图3
浮点函数异常处理方法流程图
Fig.3
Flow of floating-point function’s exception handling
算法1. ExceptionHandle:提供一个融合异常处理后的浮点函数算法,该算法分阶段对各类异常进行判断并处理.其中,COUNT_FDIV为核心运算中除法指令的个数,KeyCompute[i]为被除法指令分割开的各浮点函数子核心运算,Count统计核心运算未执行指令.其余相关变量的说明详见各子算法.
Inputs:任意浮点数构成的向量V=(p1,p2,p3,p4).
浮点函数F的定义域边界a,b;
F的参数类型pre;
F的参数个数num;
1. ExceptionCoding(E_code,E_type,E_pre,E_add); //算法2
2. if PartitionInputParameter(num,V,a,b) no INV then //算法3
3. for 1,…,COUNT_FDIV do
4. hile check_INF()≠0 and Count(KeyCompute[i]())≠do
5.R=KeyCompute[i](R);
6.nd while
7.f Count(KeyCompute())=then
8. if PartitionOutputParameter(R) no FPF then //算法5
9.eturn R; goto 3; //return result
10.else
11.=Exception(exc_num); //算法6
12.eturn E; //return exception
13.end if
14.lse
15.E=Exception(exc_num);
16. return E; //return exception
17.nd if
18.f check_DZE()≠0 do
19. E=Exception(exc_num);
20. return E; //return exception
21.nd if
21.end for
22. else
23.E=Exception(exc_num);
24.return E; //return exception
25. end if
2.1 理论依据
在对结论1进行论证前,首先关注如下的定义及运算规则,下述研究都在浮点数的背景之下进行.
定义1(机器可表示浮点数$Re $). 由0、有限数和非有限数组成的全体集合,从逻辑上可以用三元组{s,e,m}表示,即n=(-1)sxmx2e-OFFSET.其中,s表示符 位,e是指数位,OFFSET是指数偏移,m是由隐含的整数位和小数f组成的尾数.
定义2(有限数$psi $). 浮点数指数在最大值和最小值之间,且整数位恒为1的数.有限数的形式是(-1)sx(1+f)x 2e-OFFSET.其中,1是被隐含的整数位.
定义3(非有限数$bar Psi $). 除0和$psi $之外的所有浮点数,包括inf,NaN和dnp.
定义4(狭义的非有限数$hat Psi $). 除非规格化数之外的所有$hat Psi $,包括inf和NaN.
运算规则1(浮点数加/减规则). (1) 零操作数检测;(2) 比较阶码大小并对阶;(3) 尾数加/减操作;(4) 结果规格化并进行舍入处理.
运算规则2(浮点数乘/除规则). (1) 零操作数检查;(2) 阶码加/减操作;(3) 尾数乘/除操作;(4) 结果规格化并进行舍入处理.
结论1. 狭义的非有限数参与加减乘运算后的运算结果必为非有限数,即$hat Psi Theta Re bar subset bar Psi ,bar subset $表示必属于.
对于结论1的证明,将分为NaN和inf两种情况.
首先,对于NaN而言,在数学上并没有相应的数或符 与之对应,NaN通常可以理解为非法操作的一种标识,所以有$NaN odot Psi bar subset NaN,NaN odot dnpbar subset NaN,NaN odot 0bar subset NaN$(不产生DZE异常)和$NaN odot NaNbar subset NaN;$
虽然当NaN取任何值且出现$sqrt {{infty ^2} + Na{N^2}} $时,该公式的结果值为inf,但除了这种特殊情况,还是有$NaN odot infbar subset NaN$,所以无论特殊情况下结果为inf,还是一般情况下结果为NaN,都有$NaN odot infbar subset bar Psi .$
综合来看,有$NaN odot Re bar subset bar Psi $,亦有$NaNTheta Re bar subset bar Psi $.
其次,对于inf而言,先假设对于任意浮点数$x = {( – 1)^{{s_x}}} cdot {2^{{E_x}}} cdot {M_x},y = {( – 1)^{{s_y}}} cdot {2^{{E_y}}} cdot {M_y}$和z,不妨设Ex全1,My全0.即:假定y为inf,且x不为NaN.
对于加减法而言,由假设有Ex≤Ey,根据运算规则1,有$z = x + y = {( – 1)^{{s_z}}}{2^E}^{_y}({M_x}{2^{{E_x} – {E_y}}} pm {M_y}).$
令$({M_x}{2^{{E_x} – {E_y}}} pm {M_y}) = {M_z}$,则$z = {( – 1)^{{s_x} + {s_y}}}{M_z}{2^E}^{_y}$,即,浮点数z的指数依然全1.所以无论当Mz=0时z为inf,还是Mz≠0时z为NaN,都有$zbar subset hat Psi bar subset bar Psi .$
对于乘法而言,根据运算规则2,有$z = x times y = {( – 1)^{{s_z}}}{2^{{E_x} + {E_y}}}({M_x} times {M_y})$.由假设可知,无论Ex=0时Ez=Ex+Ey全1,还是Ex≠0时Ez超出可表示范围,都有$zbar subset hat Psi bar subset bar Psi $或z出现溢出.而当出现溢出时,同样可以表示为NaN,所以有$zbar subset hat Psi bar subset bar Psi .$.
综上所述,$hat Psi Theta Re bar subset bar Psi .$
结论2. 有限数与有限数进行四则运算后的结果可能为非规格化数,即$Psi odot Psi tilde subset dnp,tilde subset $表示可能属于.
对于结论2,假设y$ in
由上述假设可知:存在有限数经过某些运算后的运行结果为dnp,即$Psi odot Psi tilde subset dnp.$
通过对结论1的证明,可以充分说明分段式异常处理的有效性,即,不会出现在核心运算过程中出现异常后在计算的输出阶段不被检测的情况.换言之,核心运算中一旦出现异常,必然会将异常一步步传递下去,直到运算结束或被检测.对结论2的证明说明输出阶段对DNO异常检测的必要性,虽然在核心运算过程中不触发DNO异常(出现dnp则补指或作为0处理,具体含义见第3.2节),但在最后的计算中还是可能产生dnp(结论2),并在输出阶段触发DNO异常.
2.2 浮点异常编码
在编码的过程中,根据浮点数类型(单精度和双精度)的不同,编码也有所差异,其对应的算法如下:
算法2. ExceptionCoding:将各类异常类型编码为64位浮点数.
Inputs:E_code,错误码;
E_type,异常类型浮点函数F的定义域边界a,b;
E_pre,异常精度;
E_add,返回地址;
将64位浮点数看作由64个变量t0~t63组成的字符串,t0表示浮点数的第0位,依此类推,且初始时t0~t63均为0.
1. t0:判错误码.当错误码为EDOM时t0=0,否则t0=1.
2. t1~t6:判异常类型.当异常类型为INV时t1=1;当异常类型为DZE时t2=1;当异常类型为OVF时t3=1; 当异常类型为UNF时t4=1;当异常类型为INE时t5=1;当异常类型为DNO时t6=1.
3. t7:判返回地址.当t7=0时,表示特殊点编 的t8~t13位为返回点地址;否则,表示特殊点编 的t8~t13 位与t14~t17位作模2加操作后再作为返回点地址.
4. t8:判单双精度.当函数类型为双精度时t8=1;当函数类型为单精度时t8=0.
5. t8~t13,t14~t17:返回地址.
6. t63:异常总标识.t1~t6存在1,则t63=1,否则t63=0.
7. exc_num=t0 or (t1<<1) or (t2<<2) or … or (t63<<63);
8. return exc_num;
运行算法2实现异常类型编码,能够确保程序快速准确地判断异常类型并返回预设的异常结果.以双精度为例,经过ExceptionCoding算法后的异常编码如下:
异常标识
编码
EXC_INV_ZERO_EDOM
0x8000000000000502
EXC_INV_NANZERO_EDOM
0x8000000000004502
EXC_INV_NANINF_EDOM
0x8000000000019902
EXC_UNF_ZERO_ERANGE
0x8000000000000511
EXC_UNF_NINF_ERANGE
0x8000000000039909
EXC_DNO_DNO_ERANGE
0x8000000000000551
EXC_DZE_NINF_EDOM
0x8000000000039904
EXC_DZE_INF_EDOM
0x8000000000009504
EXC_OVF_INF_ERANGE
0x8000000000009509
EXC_ISIEEE
0xffffffffffffffff
其中,EXC_ISIEEE为特殊标识,当异常类型被设置为该标识时则直接返回,不对异常进行处理;EXC_INV_ ZERO_EDOM表示INV异常,返回值为0,且错误码为EDOM(源自于函式的参数超出范围);EXC_DNO_DNO_ ERANGE表示DNO异常,返回值为非规格化数,且错误码为ERANGE(源自于函式的结果超出范围).
2.3 输入参数检测
程序员在编程过程中,总是期望所有输入都不会引起程序的异常情况.但不幸的是:即使将输入参数小心地控制在函数定义域范围内,也可能触发一些异常(如DNO异常);幸运的是,该异常一般不会引起严重的后果.而当浮点数的特殊数参与浮点运算时,必然会引发浮点算术异常,从而使得函数运算过程中断,降低浮点函数的可靠性及性能.因此,在进入函数正常运算阶段前需要对函数的输入参数做特殊数检查:一方面将错误码为EDOM的异常解决在函数核心运算之前,另一方面确保函数核心运算的连贯性.
算法3. PartitionInputParameter:实现输入参数的预判断,确保引起EDOM类异常的输入进入相应的异常处理子程序,其他则进入核心运算子程序.
Inputes:V,输入参数向量(p1,p2,p3,p4);
a,F定义域下界;
b,F定义域上界;
num,F的参数个数1~3;
MAX,无穷数;
VCPYF($1,$2):将寄存器$1中的浮点数扩展为向量并存入寄存器$2;VINSF($1,$2,$3):将$1中的浮点数插入到向量寄存器$2中的第$3个位置;CheckINFNAN(inf,V):检测向量V中是否存在无穷数或非数;CheckInterval(a,b,V):检测向量V中的变量是否在定义域区间内.
1. VCPYF(p1,V); //V=(p1,p1,p1,p1)
2. if num=2 then
3.VINSF(p2,V,1); //V=(p1,p2,p1,p1)
4. else if num=3 then
5. VINSF(p2,V,1);VINSF(p3,V,2); //V=(p1,p2,p3,p1)
6.end if
7. end if
8. if CheckINFNAN(MAX,V)=TRUE then
9. eturn INV;
10. else if CheckInterval(a,b,V)=TRUE then
11.retuen INV;
12.else
13.return;
14.end if
15. end if
PartitionInputParameter算法的贡献并不仅仅是实现了输入参数的检测这样一个功能,而是在实现的过程中结合SIMD编程思想,确保了算法的高效性,同时也实现了多参数函数的统一处理.程序中对SIMD的运用主要体现在向量V的构建.研究中发现:浮点函数的输入参数个数从1~3不等,输入参数的不确定性,增加了标量程序的工作量,如果使用标量运算,对于3个输入的函数而言,同样的判断则需要重复进行3次;而通过向量实现,则可以同时进行判断.通过实验分析可知:与标量实现相比,向量实现平均能有10%左右的性能提升.
2.4 特定代码检测
当前,对INF异常检测的研究尤为深入,大致可分为两类.
· 一类是动态检测技术.Brumley等人[通过对可能引起INF的每一个整数操作进行检测实现的C程序检测工具RICH,该工具最大的优势在于对程序性能干扰的控制,平均只有5%的性能下降;类似的方法还有Chinchani等人[的研究.Dietz等人[针对C/C++语言并从软件工程的角度实现的检测工具IOC,并发现了众多存在的INF,该工具并没有关注于对程序性能的影响;而Chen[实现的检测方法最大的缺陷就是造成程序性能下降了50倍;
· 另一类是
声明：本站部分文章及图片源自用户投稿，如本站任何资料有侵权请您尽早请联系jinwei@zod.com.cn进行处理,非常感谢！