看完本文你将深入理解浮点数实现的原理:规格化表示、非规格化表示、+/- 0、INF、NAN、浮点数的表示范围和精度。并且也将搞清楚大数吃小数的原理以及相应的解决方案。

前几天有个同事和我讨论浮点数精度和大数吃小数的问题,正好借这个机会,把浮点数的实现原理和大数吃小数的问题写成本文,也方便后来人遇到类似问题查阅。

浮点数实现原理

这里以 IEEE-754 单精度浮点数为例,双精度浮点数也是类似的。

单精度浮点数在计算机中以 32 位存储。这 32 位被划分为符号位(Sign)、指数位(Exponent)和尾数位(Mantissa,Fraction)三个部分:

  • 符号位:表示数值的符号,即正负。占 1 bit。
  • 指数位:表示数值的指数部分。占 8 bit。
  • 尾数位:表示数值的底数部分。占 23 bit。

V = ( − 1 ) S ∗ M ∗ 2 E V = (-1)^S * M * 2^E V=(1)SM2E

在这里插入图片描述

规格化表示

尾数部分

规格化(Normal)表示首先要将小数部分做规格化处理,即表示为 1.xx * 2 ^ exp 形式。例如:

0.09375 的二进制表示为 0.00011,规格化表示为 1.1 * 2 ^ (-4)。

这里也可以看到为什么需要规格化表示:尾数位固定占 23 位,尾数位如果有许多先导零,就会占用许多有效位,使用规格化表示可以去除这些先导零,变相使有效位变多,也就提高了浮点数的表示精度

另外,规格化表示有个特点就是有效数字的首位都是 1,这个 1 也可以省去。

指数部分

指数位占 8 bit,可以表示 [0, 255] ,为了表示负数(规格化表示产生),需要减去一个偏置(bias), b i a s = 2 k − 1 − 1 bias = 2^{k-1} -1 bias=2k11,k 为指数位个数。单精度浮点的偏置为 127,则指数部分范围为 [-127, 128] ,这里去除指数位全为 0 和 全为 1 的情况(后面会介绍用于非规格化表示和特殊数),则指数部分表示范围为 [-126, 127] 。

再看上面的例子, 0.09375 ( 0.00011):

首先,规格化表示:1.1 * 2 ^ (-4),尾数部分为 1.1,省去第一个 1 ,尾数为 1。指数为 -4,加上偏置 127,等于 123。符号位为 0。在计算机中的二进制表示为:

在这里插入图片描述
用十六进制表示则为 0x3DC00000,这个大家可以使用 这个 在线工具进行验证。

在这里插入图片描述

也可以直接写代码验证下:
在这里插入图片描述

再来看一个复杂点的数值,pi = 3.1415926:

二进制表示为:11.00100100001111110110100110100010010110110000100101 。大家可以使用 这个 工具转化。

使用规格化表示,并保留小数点位 23 位:1.10010010000111111011010 * 2 ^ 1

省去前面的 1 ,则尾数位为:10010010000111111011010

指数为 1,再加上 127,等于 128,二进制表示为:10000000

符号位为 0

则在计算机中表示为:
在这里插入图片描述
对应十六进制为 0x40490FDA

最后,总结下规格化表示的公式:

M = 1.f,f 为尾数部分二进制表示。

E = exp - bias

也即规格化表示的公式:
v = ( − 1 ) S ∗ 1. f ∗ 2 e x p − b i a s v = (-1)^S * 1.f * 2^{exp-bias} v=(1)S1.f2expbias

其中 bias = 127。

对此,大家应该对浮点数的规格化表示已经非常清楚了。

非规格化表示(包括 0)

从上面介绍我们可以知道,浮点数的规格化表示是为了提高浮点数实现的精度,那为什么还需要非规格化(Denormal)表示呢?

根据浮点数的规格化表示我们知道,规格化表示能够表示的最接近 0 的数为 1.00…000 * 2 ^ (-126)。对应的二进制表示为:

在这里插入图片描述

为了能够得到更加接近 0 的数,也即扩大浮点数的表示范围,使用浮点数的非规格化表示。

比如 1.1 * 2 ^ (-130),-130 超出了规格化指数最小 -126 能够表示的范围,使用非规格化表示为 0.00011 * 2 ^ (-126),二进制表示如下:

在这里插入图片描述
因此,对于非规格化:

M = 0.f,f 为尾数部分二进制表示。

E = 1 - bias

v = ( − 1 ) S ∗ 0. f ∗ 2 1 − b i a s v = (-1)^S * 0.f * 2^{1-bias} v=(1)S0.f21bias

其中 bias = 127。

也即:
v = ( − 1 ) S ∗ 0. f ∗ 2 − 126 v = (-1)^S * 0.f * 2^{-126} v=(1)S0.f2126

非规格化表示的指数为全为 0(规格化表示指数部分不全为 0),能够表示最接近 0 的数为 0.00…01 * 2 ^ (-126)。当尾数也全为 0 ,则表示数值 0,符号位为 0,则表示 +0,符号位为 1 ,则表示 -0。

这里分析下为啥非规格化表示的 E = 1 - bias 而不是 - bias?

为了方便描述,这里以 8 bit 浮点数为例:指数位占为 4 bit,尾数位占 3 bit。则 bias = 7。

对于非规格化表示: E = 1 - bias = 1 - 7 = -6,则最大非规格化表示的最大值为 7 512 \frac{7}{512} 5127

对于规格化表示:最小规格化表示最小的数为 1.000 ∗ 2 − 6 = 8 512 1.000 * 2^{-6}=\frac{8}{512} 1.00026=5128

可见最大非规格化数和最小规格化数之间的转变会比较平滑,这得益于将非规格化时将 E 定义为 1 - bias 而不是 -bias,这样可以补偿非规格化数的有效数没有隐含的开头的 1。

特殊值(INF,NaN)

当指数位全为 1 时,尾数不全为 0,则表示一个非数,即 NaN(Not a Number)。

例如: s q r t ( − 1 ) = N a N sqrt(-1) = NaN sqrt(1)=NaN

当指数位全为 1 时,尾数全为 0,则表示一个无穷大/无穷小的数,即 INF( ∞ \infty )。符号位为 1 则为 +INF,否则为 -INF。

例如: 1 0 = + ∞ \frac{1}{0} = +\infty 01=+

表示范围

  • 单精度浮点数的表示范围:-3.40E+38 ~ +3.40E+38
  • 双精度浮点数的表示范围:1.79E+308 ~ +1.79E+308

以单精度浮点数为例,能够表示的最大数为:
( + 1.11...11 ∗ 2 127 ) 2 ≈ ( 2 128 ) 2 ≈ ( 3.4 ∗ 1 0 38 ) 10 (+1.11...11*2^{127})_2\approx(2^{128})_2\approx(3.4*10^{38})_{10} (+1.11...112127)2(2128)2(3.41038)10

精度

精度是指有效数字的最大位数,从左边第一个不为 0 的数字开始的个数。

对于单精度浮点数,由于尾数有23位,算上规格化表示小数点前隐藏的 1,共24位,24位能够表示的最大数为 2 24 − 1 = 16 , 777 , 215 2^{24}-1=16,777,215 2241=16,777,215 。最大 8 位十进制数,能够保证十进制 7 位,也即十进制的 7~8 位。

类似的,双精度浮点十进制精度为 15~16 位。

大数吃小数

看一个大数吃小数的例子:

 #include<stdio.h>
 
 int main() {
     float sum = 0.0;
     float x = 1.0;
     for (int i = 0; i < 30000000; i++) {
         sum += x;
     }
     printf("%1.10f\n", sum);
     return 0;
 }

输出结果为:16777216.0000000000

原理

要搞清楚上述问题的原因,首先要明白浮点数相加的过程:

  1. 指数位对齐。小的向大的对齐。
  2. 尾数求和
  3. 规格化

以 0.5 + 0.125 为例:

0.5 = ( − 1 ) 0 ∗ 1.0 ∗ 2 − 1 0.5=(-1)^{0}*1.0*2^{-1} 0.5=(1)01.021

0.125 = ( − 1 ) 0 ∗ 1.0 ∗ 2 − 3 0.125=(-1)^{0}*1.0*2^{-3} 0.125=(1)01.023

0.5 的指数位为 -1,0.125 的指数位为 -3。因此 0.125 的有效位需要右移,即:

0.125 = ( − 1 ) 0 ∗ 0.01 ∗ 2 − 1 0.125=(-1)^{0}*0.01*2^{-1} 0.125=(1)00.0121

然后, 0.5 + 0.125 = ( − 1 ) 0 ∗ 1.01 ∗ 2 − 1 0.5+0.125=(-1)^{0}*1.01*2^{-1} 0.5+0.125=(1)01.0121

从上面浮点数相加的过程,我们知道,由于浮点数有效位是 23 位,当较大数是较小数的 2 24 = 16777216 2^{24}=16777216 224=16777216 倍及以上时,较小数为了将指数位向较大的数对齐,有效位就会右移很多位,右移超过了 23 位,则会消失掉。

这样,上面这个例子问题的原因也就解释了。

解决方案

对于大数吃小数问题的解决思路不外乎就是使用更高精度的数据类型、尽量避免较大数和较小数直接相加,或者通过补偿的方法。

方案一:使用更高精度数据类型

使用double代替 float。

 #include<stdio.h>
 
 int main() {
     double sum = 0.0;
     double x = 1.0;
     for (int i = 0; i < 30000000; i++) {
         sum += x;
     }
     printf("%1.10f\n", sum);
     return 0;
 }

输出结果为:30000000.0000000000

方案二:分段求和

分段求和思想的本质就是避免较大数和较小数直接相加。

 #include<stdio.h>
 
 int main() {
     double sum = 0.0;
     double x = 1.0;
     for (int i = 0; i < 10000000; i++) {
         sum += x;
     }
     for (int i = 0; i < 10000000; i++) {
         sum += x;
     }
     for (int i = 0; i < 10000000; i++) {
         sum += x;
     }
     printf("%1.10f\n", sum);
     return 0;
 }

输出结果为:30000000.0000000000

方案三:Kahan 求和

Kahan 求和算法的思想:保存较大数和较小数相加过程中,较小数丢失的有效数字,然后补偿回来。Kahan 累加补偿在维基百科中有比较详细介绍,想了解的同学可以移步 这里

 #include<stdio.h>
 
 int main() {
     double sum = 0.0;
     double x = 1.0;
     float eps = 0.0;
     float t = 0.0;
     for (int i = 0; i < 30000000; i++) {
         float y = x - eps;
         t = sum + y;
         eps = t - sum -y;
         sum = t;
     }
     printf("%1.10f\n", sum);
     return 0;
 }

输出结果为:30000000.0000000000

至此,我想应该讲清楚了浮点数实现原理和大数吃小数的问题。

参考

  • [1] 深入理解计算机系统
  • [2] https://en.wikipedia.org/wiki/Kahan_summation_algorithm
Logo

开放原子开发者工作坊旨在鼓励更多人参与开源活动,与志同道合的开发者们相互交流开发经验、分享开发心得、获取前沿技术趋势。工作坊有多种形式的开发者活动,如meetup、训练营等,主打技术交流,干货满满,真诚地邀请各位开发者共同参与!

更多推荐