1. 非规则浮点数(subnormal)
根据IEEE754(http://zh.wikipedia.org/wiki/IEEE_754),浮点数由符号位、指数、尾数组成,对8字节的double 类型数据,指数位有11位,可以表示2的阶数范围 -1022~+1023,对应的指数部分的值是1~2046(加上1023的偏移)。
double类型的尾数位有52位,尾数默认的最高有效位为1,存储时省略,所以规则的浮点数范围是 +-2^-1022 ~ +- (2-2^-52)*2^1023,最小的规则正数大约是4.5E-307。除去0以外,指数部分取值0(尾数不为0)表示一类特殊的数——subnormal,取值小于2^-1022。
2. CPU对subnormal 数的处理
现代CPU对subnormal一般是硬件处理,ALU只能直接运算normal 数,对subnormal 触发处理器异常,然后处理。以Intel CPU为例,在Intel 的手册中,subnormal 写作denormal,处理涉及两个方面:
(1) normal 数的运算结果是 denormal:CPU设置 numeric underflow exception(#U),通过称作“gradual underlow”的运算过程,对尾数逐位向右移位,直到指数增大到正常范围。
(2) denormal 是源操作数:CPU设置 denormal operand exception(#D),这个异常在计算指令真正执行之前发生。处理过程与寄存器MXCSR中指定的相关模式有关,默认是与IEEE754兼容,但是也可以通过设置“denormal-are-zeros”模式(奔腾4和Xeon),将denormal操作数置0提升运算速度。
以上的处理器异常默认是由CPU自己处理(即masked),但是也可以通过设置,将这些异常暴露给用户处理。
3. subnormal 运算的性能
这里可以用一个microbenchmark来测试一下,不同的操作数,乘法运算效率相差有多少:
04 | static void mul_init( double x ) { |
06 | __asm__ __volatile__ ( |
07 | 'movsd %1, %0' : '=Yz' (dummy): 'x' (x) |
09 | __asm__ __volatile__ ( |
16 | __asm__ __volatile__ ( |
17 | 'movsd %1, %0' : '=Yz' (x): 'x' (dummy) |
19 | __asm__ __volatile__ ( |
26 | static void mul_latency_pack4() { |
27 | __asm__ __volatile__ ( |
35 | static void mul_throughput_pack4() { |
36 | __asm__ __volatile__ ( |
44 | static void measure_latency( int repeats ) |
48 | REPEAT_256( mul_latency_pack4(); ) |
53 | static void measure_throughput( int repeats ) |
57 | REPEAT_256( mul_throughput_pack4(); ) |
62 | int main( int argc, char *argv[]) |
69 | TIMING( measure_latency(repeats/10), cost ); |
70 | printf ( 'Latency (normal): %.3e ns\n' , cost * 1e9 * 10 / ( double )repeats / 1024 ); |
73 | TIMING( measure_throughput(repeats), cost ); |
74 | printf ( 'Throughput (normal): %.3e ops\n' , ( double )repeats * 1024 / cost ); |
76 | double fsubnormal = 1e-310; |
77 | mul_init( fsubnormal ); |
80 | TIMING( measure_latency(repeats/10), cost ); |
81 | printf ( 'Latency (subnormal): %.3e ns\n' , cost * 1e9 * 10 / ( double )repeats / 1024 ); |
84 | TIMING( measure_throughput(repeats), cost ); |
85 | printf ( 'Throughput (subnormal): %.3e ops\n' , ( double )repeats * 1024 / cost ); |
编译参数 -O3 -finline-functions,在我的2.4GHz主频的机器上,一个结果是:
1 | Latency (normal): 2.090e+00 ns |
2 | Throughput (normal): 1.916e+09 ops |
3 | Latency (subnormal): 7.111e+01 ns |
4 | Throughput (subnormal): 1.871e+07 ops |
相对于规则的浮点数,subnormal 运算的延迟增大到34倍,吞吐量减小100倍。
4. 解决方法
科学计算里,subnormal 有时很难避免。除了运算变慢之外,由于尾数的有效位减少,精度也降低,后续的运算有产生无穷大inf的风险等等。处理的难度:
(1) 无法确定哪一步操作将产生subnormal 时,对所有中间结果检查开销很大。
(2) subnormal 仍然是有意义的数,直接置0是不是会产生错误的结果(比如除0错误)。
对于inf 和nan 类的不正常的数,也同样有上面类似的问题。
理想的情况是,从算法上做一些处理,比如对产生subnormal 的源操作数做偏移或者截断,消除源头。总之还是比较棘手的。