深入剖析浮点存储及其运算规则

 

对于浮点类型数据，首先我们需要明白的一点是：浮点数和整型数的编码方式是很不一样的，IEEE浮点标准采用V = (-1)^s×M×2^E的形式来表示一个数，其中符号s决定是负数(s=1)还是正数(s=0)，由1位符号位表示。有效数M是一个二进制小数，它的范围在1~2-ε之间（当指数域E既不全为0也不全为1，即浮点数为规格化值时。ε为有效数M的精度误差，比如当有效数为23位时，ε为2^-24），或者在0~1-ε之间（当指数域全为0，即浮点数为非规格化值时），由23位或52位的小数域表示。指数E是2的幂，可正可负，它的作用是对浮点数加权，由8位或者11位的指数域表示。

 

下面我们来详细地剖析IEEE浮点数据表示，相信在认真研读完以下内容之后，你对浮点数的存储将会有非常清晰的认识。不用担心文中会充满数学家才会考虑的算式和公式，虽然大多数人认为IEEE浮点格式晦涩难懂，但理解了其小而一致的定义原则之后，相信你能感觉到它的优雅和温顺，以下的叙述将尽量用浅显易懂的语言和例子让大家心情舒畅地了解浮点存储的内幕。

 

以32位浮点数为例，其存储器的内部情况是这样的：

 

   符号域               指数域                                                     小数域

 

S(1位)

E(8位)

M(23位)

 

从上述公式V = (-1)^s×M×2^E可以计算出一个浮点数具体的数值，要理解这三个数据域是如何被解释的，我们需要知道：浮点数的编码根据指数域的不同取值表示被分成三种情形：

1、规格化值。

当指数域不全为0且不全为1的时候即为这种情形，这是最常见的情况。这时有效数域被解释为小数值f，且 0≤f＜1，其二进制表示为0.f_n-1f_n-2…f₁f₀，而有效数M被解释为为M=1+f。同时指数域被解释为偏置形式的有符号数，就指数E被定义为E=e-Bias，其中e就是指数域的二进制表示，Bias是一个等于2^k-1-1的偏置值(k为指数域位数，对于32位浮点数而言是8，对于64位浮点数而言是11)。我们可以调整指数域E来使得有效数M的范围在1≤M＜2之间，也就可以始终使得M的第一位是1，从而也没有必要在有效数域中显式地表示它了。

 

我们来做一个透视浮点数存储的练习，以此来更好地理解以上内容。例如浮点数12345.0，其二进制表示为(1.1000000111001)₂×2¹³，显然符号位S应该为0，而指数域部分E应该根据公式E= e-Bias计算，e就是我们想要的内部存储表示，Bias在32位单精度浮点中的值为127，因此e=E+Bias=13+127=140，用二进制表示即为(10001100)₂，而小数域就是在其二进制表示的基础上减去最高位的1即可，即0.1000000111001，补满23位小数位，即可得到0.10000001110010000000000，存储时略去前面的的整数部分和小数点，因此浮点数12345.0的IEEE浮点表示为：

 

0

1000 1100

100 0000  1110 0100  0000 0000

 

2、非规格化值。

当指数域全为0时属于这种情形。非规格化编码用于表示非常接近0.0的数值以及0本身（因为规格化编码时有效数M始终大于等于1），我们会看到，由于在非规格化编码时指数域E被解释为一个定值：E = 1-Bias = -126（而不是规格化时E = e-Bias），使得规格化数值和非规格化数值之间实现了平滑过渡。另外，此时有效数M解释为M=f，对比上面所讲的规格化编码，此时的有效数M就是小数域的值，不包含开头的1。

 

举个例子，编码为如下情形的浮点数就是一个非规格化的样本：

0

0000 0000

000 0000  0000 0000  0000 0001

 

这个数值的大小，就应该被解释为(-1)⁰×(0.00000000000000000000001)₂×2^1-Bias，即

2^-23×2^-126 = 2^-149，转成十进制表示大约等于1.4×10^-45，实际上这就是单精度浮点数所能表达的最小正数了。

以此类推，规格化值和非规格化值所能表达的非负数值范围如下所示：

 

 

指数域

小数域

32位单精度浮点数

值

十进制

0

0000 0000

000 0000 0000 0000 0000 0000

0

0.0

最小非规格化数

最大非规格化数

最小规格化数

最大规格化数

0000 0000

0000 0000

0000 0001

1111 1110

000 0000 0000 0000 0000 0001

111 1111 1111 1111 1111 1111

000 0000 0000 0000 0000 0000

111 1111 1111 1111 1111 1111

2^-23×2^-126

(1-ε)×2^-126

1×2^-126

(2-ε)×2¹²⁷

≈1.4×10^-45

≈1.2×10^-38

≈1.2×10^-38

≈3.4×10³⁸

1

0111 1111

000 0000 0000 0000 0000 0000

1×2⁰

1.0

规格化和非规格化值非负数值范围

 

从上表中可以看出，由于在非规格化中将指数域数值E定义为E = 1-Bias，实现了其最大值与规格化的最小值平滑的过渡（最大的非规格化数为(1-ε)×2^-126，只比最小的规格化数2^-126小一点点，ε为2^-24）。

 

2、特殊数值。

当指数域全为1时属于这种情形。此时，如果小数域全为0且符号域S=0，则表示正无穷+∞，如果小数域全为0且符号域S=1，则表示负无穷-∞。如果小数域不全为0时，浮点数将被解释为NaN，即不是一个数(Not a Number)。比如计算负数平方根或者处理未初始化数据时。

 

以下是理清各种数据之间关系的总结：

1.  浮点数值V = (-1)^s×M×2^E。

2.  在32位和64位浮点数中，符号域s均为1位，小数域位数n分别为23位和52位，指数域位数k分别为8位和11位。

3.  对于规格化编码，有效数M = 1+f，指数E = e-Bias（e即为k位的指数域二进制数据，对于32位浮点数而言e的范围是1~254，此时E的范围是-126~127）。

4.  对于非规格化编码，有效数M = f，指数E = 1-Bias（这是一个常量）。

5.  Bias为偏置值，Bias = 2^k-1-1，k即为指数域位数，在32位和64位浮点数中k分别为8和11。

6.  f为小数域的二进制表示值，即n位的小数域f_n-1f_n-2…f₁f₀将被解释为f=(0.f_n-1f_n-2…f₁f₀)₂。

 

 

有了以上的背景知识之后，我们就可以更从容地分析浮点运算了。毕竟我们不是数学家，学习浮点数不是为了科学研究，更多地是从实用主义的角度出发，是为了要写出更好更可靠的代码。

 

首先要说明的是，浮点运算是不遵循结合性的，也就是说(a+b)-c和a+(b-c)可能会得出不一样的结果（比如(1.23+4.56e20)-4.56e20 = 0，但1.23+(4.56e20-4.56e20) = 1.23）。这是由浮点数的表示方法限制的，浮点数的范围和精度有限，因此它只能近似地表示实数运算。在两个浮点数进行运算的时候，它们之间会产生一种称之为“舍入”的行为，要理解清楚这种游戏规则，请仔细研读以下叙述。

 如果现在要执行f1 + f2，则有可能会发生舍入操作，具体情况如下：

假设有两个浮点数f1和f2，其IEEE存储格式分别为如下：

s

e_k-1e_k-2e_k-3…e₁e₀

f_n-1f_n-2f_n-3……f₁f₀

 

s’

e’_k-1e’_k-2e’_k-3…e’₁e’₀

f’_{n-1f’n-2f’n-3}……f’₁f’₀

对于单精度浮点数而言k=8，n=23，对于双精度浮点数而言k=11，n=52。这个对于我们讨论舍入问题不是关键，权且就把f当成是32位的单精度浮点数即可。

 

1、 取两数中指数域较大的一个e，再取小数域为最小值0000…0001。

2、 令δ= M×2^E，即(0.0000…0001)₂×2^e-127，注意：这里M取的是非规格化值。

3、 所有小于δ/2的数值都会被舍入。

 

 

下面的代码印证了上面的推论。

// example.c#include <stdio.h> int main(void){        float  x = 0.1;        float  a = x + 0.37252e-8;        float  b = x + 0.37253e-8;         if( x == a )                printf(“x == a/n”);        if( x == b )                printf(“x == b/n”);         return 0;}

程序的运行结果是打印出了x==a，但是不会打印x==b。也就是说系统辨识不出来比0.372529e-8还小的值，区分不出来x与a的差别。其内幕如下：

单精度数据x=0.1的存储细节如下：

0

011 1101 1

100 1100 1100 1100 1100 1101

因此δ/2 = (M×2^E) /2 = (2^-23×2^123-127)/2 = 2^-28≈ 0.372529e-8。变量a与b都分别在x的基础上加了一个常量，但是由于加在a上的增量小于δ/2，在此精度范围内该增量无法被辨识。由此可见，从数学角度上看明明是三个不相等的实数，但是由于计算机本身特性的限制会导致程序运行结果跟预料的不符，这就要求我们必须对浮点数运算的内部实现非常了解，不能想当然地写出似是而非的代码。

同样地，我们可以自己再做一个练习，比如有一个单精度浮点数f = 10e18，用以下代码查看其二进制存储细节：

float f = 10e18; FILE  *fp  =  fopen(“abc”,  “w+”); /* 产生一个叫abc的文件来存储f */fwrite(&f,  sizeof(f),  1,  fp); /* 将f以二进制形式（而非字符形式）写入该文件 */

运行命令 od  -x  abc查看其内部细节

 

然后计算出δ/2≈ 0.549755e12，也就是说对于f而言精度小于0.549755e12的数据都将会被舍入，所有跟f的差值小于此精度的数据都将会被认为等于f。

 

从上述叙述中我们得出一个重要的结论：对于每个不同的浮点数，都有相应的最小可辨识精度（即δ/2），此最小可辨识精度随着该浮点数的数值变化而变化，具体究竟是多少要像上面那样分析该浮点数的二进制存储内部细节，找到其指数域之后才能确定，我们根据这个最小可辨识精度才能明确判定代码中所有对此浮点数的运算是否有效，否则可能会由于舍入的问题存在而在逻辑上存在歧义。

 

网上有文章简单地用一个固定的最小精度EPISON来判断两个浮点数是否相等的方法是不严谨的。我们以后在程序中需要进行浮点比较的时候，最可靠的方法是先计算出其相应的最小可辨识精度δ/2，在此基础上再进行操作和判断。