浮点数
最近有个朋友碰到一个奇怪的问题,两个浮点数相加得到的结果并非数学意义上的正确结果,对于这一点也许大家都清楚计算机并不能精确地表示浮点数,这就是所谓的浮点数的精度问题。但是至于为什么计算机不能精确地表示浮点数?两个两位小数加法运算的结果为什么得到的结果会有很多位小数?对于这两个问题大家通常都只是简单把它归结为浮点数的精度问题导致的。如果要彻底搞清楚为何会如此,我们就需要真正理解计算机是如何表示和存储浮点数的,以及计算机中浮点数的运算规则,这也是本文所要探讨的主题。
首先看一个简单的示例,这个示例来自于StackOverflow上的一个帖子:
<?php
$a = 35;
$b = -34.99;
echo ($a + $b); //0.009999999999998
在我的机器上运行得到的结果为0.009999999999998,我的机器是64位机器,Mac OS系统,PHP也是64位的。据说有些机器和系统中可以得到精确的结果,这可能是有些系统对浮点数的运算做了一些特殊处理,例如Windows,这不属于本文的讨论范围。现在绝大部分机器的浮点数的表示和运算都是遵从IEEE745规范,我们的讨论也是基于这个规范的,对于希望深入了解浮点数的同学可以自行搜索这个规范。
科学记数法
科学记数法应该是初中代数中的内容,它的概念很简单,下面是一个来自百度百科的定义:
科学记数法是一种数学专用术语。将一个数表示成 a×10的n次幂的形式,其中1≤|a|<10,n为整数,这种记数方法叫科学记数法。例如920000可以表示为9.2*105,读作9.2乘10的5次方。
这个定义针对的是十进制的实数,根据这个定义,所有的十进制的实数都可以使用科学记数法表示。同样肯定也存在二进制数的科学记数法,例如1011 1101可以表示为:
1011 1101 = 1.011 1101 * 2^7
从这个示例中可以看出科学记数法中包括三个部分:
- 有效数,在这个示例中是1.011 1101
- 基数,二进制的基数是2,十进制的基数是10
- 幂,在这个示例中是7
有效数可以是正数,也可以是负数。正确的科学记数法对有效数是有限制的,对于二进制而言,它必须大于等于1,并且小于2。因为二进制表示中只有0和1,所以这个限制其实也是规定了有效数的整数位只能为1,后面讨论计算机中的表示的时会再次谈到这点。
上面的示例中所表示的二进制数实际上是一个整数,而我们要讨论的是浮点数,在数学上也许浮点数可以定义为存在小数位的数,但是在计算机中却并非如此。因为计算机中整数的表示和浮点数的表示是不同的,就以上面提到的php代码为例,$a在定义的时候是一个整数,但它在跟$b相加的时候会先转换为浮点数,这些都是编译器和运行环境负责处理的。
浮点数的存储
学过C语言的都知道C中有两种表示浮点数的类型float和double,我们把float称为单精度浮点数,而double则是双精度浮点数,它们的最主要差别是float类型的变量所占的内存空间是4个字节,每个字节包含8个二进制位,所以4个字节会有32个二进制位;而double则需要8个字节的内存空间,64个二进制位。这个差别就是所谓的精度的差别,精度越高它可以表示的浮点数的范围越大,而且浮点数的位数会越多,关于这一点后面会进一步阐述。
很多机器中还会有一个双精度扩展(double extended)。在C中它的类型是long double,对于32位机器,这个类型跟double类型没有差别,也是8个字节。但是在64位机器上,它会占用10个字节,也就是80位(实际上是79位,有一位是冗余的),实际上在以4个字节对齐的机器上它所占用的内存空间是12个字节。双精度扩展不属于本文的讨论重点,我们后面不会在谈到它。
单精度浮点数
IEEE规定单精度浮点数由三个部分组成:
- 小数部分,它包含23个二进制位,我们把它命名为f,它实际上是科学记数法中有效数的小数位
- 幂,它包含8个二进制位,我们把它命名为e
- 符号位,它只有一个二进制位,所以它的值只可能是0或1,我们把它命名为s
对于任一浮点数N,我们可以用下面的公式来定义N:
N = (-1)^s * 1.f * 2^e
这里要重点注意一下1.f,我们上面说过二进制的科学记数法中有效数的整数位总是为1,所以在计算机中,这一位是作为保留位,它不会被存储。23位的f,8位的e和1位的s,它们加起来就是32位,刚好存储在4个字节的内存空间中,下图就是这个4字节的内存空间的示意图:
在这个图中,0位表示最低位(least significant),0到22位存储的是f,23到30位存储的是e,31位存储的是s。对于幂的存储,IEEE会有些特殊处理,它会先把真正的幂加上127再存储在23到30位的8个二进制位中。如果我们把这8个二进制位表示的数定义为e,那么真正的幂是e-127,所以如果我们是从内存中读取s、e、f,那么上面表示N的公式将修改为:
N = (-1)^s * 1.f * 2^(e-127)
实际上IEEE除了定义了正常的浮点数的表示外,还定义了一些其他特殊的数的表示,对于正常数e的范围为:0 < e < 255,这些特殊数包括非正常数(低于正常数)、0(正数0和负数0)、正无穷、负无穷以及NaN(Not-a-Number),后面会这些数作进一步的说明。下面我们先探讨怎么把一个十进制的浮点数转换为一个二进制的单精度浮点数,从而可以把它存储在内存中。
还有一点需要强调下,23位所存储的小数部分,加上整数1,也就是1.f总共有24位,这24位也被定义为精度位,32位的单精度浮点数有时候也被称为提供24位精度的浮点数。
十进制浮点数转换为二进制单精度浮点数
这是一个很有意思的问题,我个人觉得这个问题并不适合用文字来讲解,我在youtube上面找到了两个视频,这个和这个,当然你需要能够翻墙才能观看。这两个视频所讲解的转换方法都不一样,第一种方法介绍了一种转换公式,确切地说是一系列公式,通过这些公式计算出幂和一个整数,这个整数转换为二进制格式就是一个24位的数,最高位是1,后面的23位就是f,这个计算使用了对数运算。第二个视频则是通过一种逼近算法来计算出幂以及f,这个方法更好理解,但有些繁琐。这里我会使用第二个视频介绍的方法来转换一个十进制浮点数。
我们要转换的十进制浮点数为936.35。我们先把它转换为1.xxxxxx * 2^m,这只需要不断除以2,直到余数处于[1,2)的范围内。整个计算过程如下:
936.35 / 2 = 468.175 2^1
468.175 / 2 = 234.0875 2^2
234.0875 / 2 = 117.04375 2^3
117.04375 / 2 = 58.521875 2^4
58.521875 / 2 = 29.2609375 2^5
29.2609375 / 2 = 14.63046875 2^6
14.63046875 / 2 = 7.315234375 2^7
7.315234375 / 2 = 3.6576171875 2^8
3.6576171875 / 2 = 1.82880859375 2^9
很显然1.82880859375转换为二进制数的结果肯定还是一个[1,2)区间内的数,所以9就是我们所要求出的幂值。下一步就是把0.82880859375转换为二进制形式,这个等会再说,先看下把二进制浮点数怎么转换为十进制浮点数。假设一个n位的二进制小数,包含i位小数,整数部分有n-i位。我们以下面的形式表示这个数,其中bn表示第n位的二进制位(0或者1):
bn-1 … bi.bi-1 …. b0 = bn * 2^(n-i-1) + bn-1 * 2^(n-i-2) … + bi * 2^0 + bi-1 * 2^(-1) + bi-1 * 2^(-2) + … + b0 * 2^(-i)
根据这个公式,假设有一个二进制浮点数:101.01011,对于这个数n=8,i=5,它的十进制表示的值为:
101.01011 = 1 * 2^(8 - 5 -1) + 0 * 2^(8–5-2) + 1 * 2^ (8 - 5 - 3) + 0 * 2^(-1) + 1 * 2^(-2) + 0 * 2^(-3) + 1 * 2^(-4) + 1 * 2^(-5)
= 1 * 2^2 + 0 * 2^1 + 1 * 2^0 + 0 * 2^(-1) + 1 * 2^(-2) + 0 * 2^(-3) + 1 * 2^(-4) + 1 * 2^(-5)
= 4 + 0 + 1 + 0 + 1/4 + 0 + 1/16 + 1/32
= 5 + 11/32
= 5.34375
根据这个公式,我们只需要对十进制浮点数的小数部分进行乘2运算,得到的结果的整数位的值就是当前二进制位的值,然后再取这个结果的小数部分再次进行乘2运算,得到的结果的整数位就是下一个二进制位的值,依次进行下去。我们先看下5.34375中的小数部分0.34375的二进制表示的转换过程:
0.34375 * 2 0.6875 0
0.6875 * 2 1.375 1
0.375 * 2 0.75 0
0.75 * 2 1.5 1
0.5 * 2 1 1
0 0 0
最终的结果为01011,跟上面的101.01011的小数位部分一致,上面的运算过程在运算结果为0的时候终止。我们再看下936.35的情况,对于936.35的情况我们要转换的十进制小数位是0.82880859375,它的计算过程为:
0.82880859375 * 2 1.6576171875 1
0.6576171875 * 2 1.315234375 1
0.315234375 * 2 0.63046875 0
0.63046875 * 2 1.2609375 1
0.2609375 * 2 0.521875 0
0.521875 * 2 1.04375 1
0.04375 * 2 0.0875 0
0.0875 * 2 0.175 0
0.175 * 2 0.35 0
0.35 * 2 0.7 0
0.7 * 2 1.4 1
0.4 * 2 0.8 0
0.8 * 2 1.6 1
0.6 * 2 1.2 1
0.2 * 2 0.4 0
0.4 * 2 0.8 0
0.8 * 2 1.6 1
0.6 * 2 1.2 1
0.2 * 2 0.4 0
0.4 * 2 0.8 0
0.8 * 2 1.6 1 //0110会无限循环下去
0.6 * 2 1.2 1
0.2 * 2 0.4 0
最终没有出现为0的情况,但是会出现一种无限循环的情况。因为在计算机中能够表现的位数是有限的,对于我们现在讨论的单精度浮点数,它只能保存23位的小数,所以936.35的单精读浮点数表示的f部分为:
110 1010 0001 0110 0110 0110
现在我们也可以很直观地看出计算机不能精确地表示这个数,这就是所谓的精度问题存在的原因。
我们现在搞清楚了怎么把一个十进制浮点数转换为二进制浮点数。对于936.35,经过转换得到它的幂为9,所以e=9+127=136,它的二进制表示为1000 1000,而f=110 1010 0001 0110 0110 0110,它是一个正数,所以s=0,所以用科学记数法表示为:
936.35 = 1.110 1010 0001 0110 0110 0110 * 2^9
在计算机的存储中,0到22位会存储110 1010 0001 0110 0110 0110,23到30位会存储1000 1000,而符号位则会存储0(正数为0,负数为1),整个4字节的内存块看起来就是下面这个样子:
0 1000 1000 110 1010 0001 0110 0110 0110
特殊数
上面我们谈到,在单精度浮点数的二进制表示中,当它的幂(加了127后的值,也就是内存中保存的值,这是e的值)处于(0,255)范围内时,IEEE规定它表示的是正常数(normal number),IEEE除了规定了正常数外,还有一些其他的数:非正常数(subnormal number,也叫denormalized number,比正常数更小的数),0(无符号和有符号),正无穷和负无穷,以及NaN,它们的定义如下表(u表示可以为任意值):
从这个表格可以看到,当e<255,它可以表示正常数,非正常数和0。注意存储e的8位二进制数所表示是一个无符号数,它的范围是[0,255]。另外所谓非常正常数其实是极小的数。当0<e<255时,计算机才会使用一个为1的保留位作为有效数的整数部分。这点实际上这个也很好理解,结合我们上面说的科学记数法的定义,有效的科学记数法对有效数的整数位是有限制的,满足这个限定的数才是正常数,而不满足这个限定的就是非正常数(或者是0),表现在计算机中就是保留位是为1还是为0。
下面的表格列出了几个单精度浮点数的示例,它现实出了每个数的位模式(以16进制表示)和对应的十进制数。其中最大的正值(大于0的)正常数是IEEE规定的最大的单精度有穷浮点数。而最小的正值(大于0的)非正常数是IEEE规定的最小的正值单精度浮点数,最小的正值正常数通常被认为是发生向下溢出的界限(underflow threshold)。另外还有一点需要指出的是,最小正值正常数跟最大正值非正常数非常近似。
还有一个需要特别指出的是,这个示例中的NaN并不是唯一的,只需要满足我们上面的表格中NaN的定义的数都可以表示NaN,这里只是选了一个作为示例。
双精度浮点数
双精度浮点数也是由三部分组成,只是每个组成部分所占的二进制位数稍有不同,小数部分有52位,幂部分有11位,符号位还是1位,总共占64位,8个字节,这三个部分的示意图如下图所示:
下表是双精度浮点数各种特殊数的定义:
这个跟单精度的情况类似,正常数幂的范围是(0,2047),对于双精度浮点数,保存在内存中的幂值等于真正的幂加上1023,另外双精度浮点数的精度位一共有53位。下面的表格是一些双精度浮点数的示例:
也可以使用上面介绍的方法把一个十进制的浮点数转换位双精度的二进制浮点数,下一步我们以本文开头的php代码中的计算为例讲解下浮点数的运算。
浮点数的运算
我们已经分别讨论了浮点数的存储和表示,以及怎么把十进制的浮点数转换为二进制浮点数从而可以在计算机中存储,所有这些工作都是为了能够在计算机中计算浮点数,这就是我们这一节要讨论的重点。鉴于篇幅关系,我们只讨论浮点数的加法(减法)运算,而乘法和除法运算可以转换为加法和减法运算来实现,并且计算机也是这么做的。我们使用双精度浮点数来表示和计算。目前64位计算机已经非常普及了,而且php中根本就不存在单精度浮点数,虽然使用64位会有点麻烦,但显然这会更接近现实情况。
我们要计算的数就是本文开头中php代码中的运算35+(-34.99),这个运算可以转换为35-34.99,所以我们先看看35和34.99的浮点数表示,这里我直接给出结果,略去转换过程,这个大家可以自行进行:
35 = (-1)^0 * 1.00011 * 2^5
34.99 = (-1)^0 * 1.000 1011 1111 0101 1100 0010 1000 1111 0101 1100 0010 1000 1111 1 * 2^5
这两个数的幂相同,幂在内存中的表示为5+1023=100 0000 0100,符号位都为0,它们在内存中的表示分别为:
35: [0] [100 0000 0100] [000 1100 0000 0000 0000 0000 0000 0000 0000 0000 0000 0000 0000 0]
34.99: [0] [100 0000 0100] [000 1011 1111 0101 1100 0010 1000 1111 0101 1100 0010 1000 1111 1]
-: [0] [100 0000 0100] [000 0000 0000 1010 0011 1101 0111 0000 1010 0011 1101 0111 0000 1]
上面也给出了相减的结果,它的幂不会变,但是有效数的整数位现在为0,如果要表示为有效的科学记数法的格式,需要做一些调整,调整后的结果为:
35 - 34.99 = (-1)^0 * 1.010 0011 1101 0111 0000 1010 0011 1101 0111 0000 1 * 2^-7
这个数在内存中的表示为:
[0] [011 1111 1000] [010 0011 1101 0111 0000 1010 0011 1101 0111 0000 1000 0000 0000 0]
把这个数转换为十进制表示得到的值为:
9.99999999999801048033987171948E-3
这跟我们程序运行的结果0.009999999999998刚好一致。CASE CLOSED:)
一些其他相关的东西
1. IEEE定义了5种浮点数异常,它们分别是:无效操作(invalid operation) 、除0(division by zero)、溢出(overflow)、向下溢出(underflow)、不准确(inexact),IEEE也规定了这些异常发生时的条件码,据我所知,CPU会处理“除0”这个异常,并且用户可以定义相关操作来处理这些异常,这篇文章对此有详细的说明。
2. 有一篇非常经典的关于浮点数理论方面的论文:What Every Computer Scientist Should Know About Floating-Point Arithmetic,基本上大多数关于浮点数的文章都会引用这篇论文。我第一次知道这篇论文是在一篇关于内存的论文(What Every Programmer Should Know About Memory)中看到的。你可以发现这两篇论文的标题很相似,后者是模仿前者的,并且后者在文中说道如果你没有看过这篇关于浮点数的论文,那就别把手放到键盘上,言外之意就是没看过那篇论文就别说你会写程序,哈哈,有些夸张。我只能说也没读完,太艰深了。
3. 浮点数的运算比起整数运算更耗性能,通过这篇文章也可以看出个一二来,所以我们在写程序的时候要尽可能避免使用浮点数。操作系统的内核都没有使用浮点数,一般用户进程进行系统调用的时候,操作系统是不会保存浮点数相关的寄存器的,因为它不会用到它们。
4. 在有专门的浮点数寄存器之前,浮点数是放到内存中的,对于浮点数寄存器相关的懂东西我不太了解,如果想了解跟浮点数相关的硬件的同学可以参考这篇文章,这个是x86的,因为这几篇文章重点是谈论SPARC的,你也可以从这篇文章中了解SPARC中跟浮点数相关的硬件。其实如果你随便写一个简单的包含浮点数变量的C程序,把它编译成汇编代码,你会发现这些变量是保存在内存中的。
5. 浮点数的操作指令跟整型数是完全不同的,有兴趣的同学可以去了解下x86系列的浮点数操作指令。
6. 这篇文章重点参考了这篇文章,也是上面列出的几篇文章所属的系列的一篇,这篇文章还讲解了双精度扩展以及向下溢出,希望更进一步了解这些知识的同学可以好好读一下这篇文章。