浮点数 Kahan求和

缘起

在计算机组成与体系结构中，浮点数是相当重要的存在. 很多超算性能测试中——每秒浮点运算次数是一个相当重要的指标. 浮点运算性能可以直观地反映一个cpu的计算能力，注意是“计算能力”，可是学过编程的人都知道，占代码量80%的是由if ，while， for 等等构成的分支语句，这些语句对cpu的浮点运算要求不高，可以说没什么要求，但要求有大量的分支预测机制，以加快速度。真正对浮点要求高的是视频压缩，场景的渲染，光线散射的计算等等。那么计算机中对浮点数的存储是怎么样的呢? 以及所谓的精度丢失是怎么回事? 如何解决? 我们来谈谈这个话题. 本文汲取了【1】~ 【5】中的营养.

分析

首先为什么会有浮点数，或者说怎么理解浮点这个词儿~ 提到浮点数，就不能不提到定点数.

早期的计算机使用定点数来表示实数，由于定点数的小数点位置固定，而计算机字长有限，定点数无法表示很大和很小的实数. 所以浮点数应运而生. 什么是字长? 我的天，你的计算机基础需要补补啦~ 先别忙着CRUD啦~

1
2
3
4
5
6
7
8
9
10
11
12
计算机字长就是计算机中CPU在一次操作中能处理的单位字的长度，即运算器能够并行处理和存储器每次读写操作时能
包含的二进制码的位数。
字长是直接用二进制代码指令表达的计算机语言，指令是用0和1组成的一串代码，
它们有一定的位数，并分成若干字长段，各段的编码表示不同的含义，
例如某台计算机字长为16位，即有16个二进制数合成一条指令或其它信息。
16个0和1可组成各种排列组合，通过线路变成电信号，让计算机执行各种不同的操作.

说的更直白一点

计算机进行数据处理时，一次存取、加工和传送的数据长度称为字（word）。一个字通常由一个或多个（一般是字节的整数位）字节构成。例如286微机的字由2个字节组成，它的字长为16；486微机的字由4个字节组成，它的字长为32位机。
计算机的字长决定了其CPU一次操作处理实际位数的多少，由此可见计算机的字长越大，其性能越优越。
现在的64位机的CPU字长为： 64 位

而浮点数在计算机中是如何存储的呢? 这就不得不提到现在最为流行的IEEE754标准（现在的浮点数都是基于此标准进行存储的）

对于规约化浮点数，IEEE754规定它的存储模型是

至于什么是规约化浮点数，什么是非规约化浮点数，以及为什么要有非规约化浮点数，且听我娓娓道来.

为什么叫规约化? 即符合一定的标准——规约的英文就是regular. 这里的regular体现就体现在有效位数M上——它必须要在[1,2)范围内. 而非规约化浮点数就是打破了这一点. 后面会详细哔哔的.

不过，如果真的要搞清楚浮点数是如何存储的，首先你得知道给你一个小数，如何将它转换为二进制. 举个栗子吧

栗子1: 将176.0625表示为二进制数

首先将176转换为二进制——10110000

​ 176/2=88…0（余数）
88/2 =44…0
44/2 =22…0
22/2 =11…0
11/2 =5 …1
5/2 =2 …1
2/2 =1 …0
1/2 =0 …1
将余数按从下往上的顺序书写就是：10110000

然后将小数部分0.0625转换为二进制, 转换方法和整数部分不断除以2取余数部分相反, 而是不断的乘以2并取整数部分

0.0625*2=0.125…0（整数）
0.125*2 =0.25 …0
0.25*2 =0.5 …0
0.5*2 =1.0 …1
将整数按从上往下的顺序书写就是：0001

所以176.0625转换为二进制就是 10110000.0001

栗子2: 将9.0 转换为二进制.

用上面的方法就是 1001.0

栗子3: 将-9.0转换为二进制.

用上面的方法并遵守IEEE754，就是 -1001.0

好了，我们知道了一个小数该怎么转换为二进制小数. 那么该怎么将这个二进制小数存储在计算机呢? 还是IEEE754规定

  1. IEEE 754规定，对于32位的浮点数(单精度浮点)，最高的一位是符号位s，接着的8位是指数E，剩下的23位为有效数字M
  2. 对于64位的浮点数（双精度浮点）来说，最高的一位仍为符号位s，接着的11位是指数E，剩下的52位为有效数字M

下面具体讨论一下 S，E，M的存储

首先S的存储没什么好说的，1就表示该浮点数是负数，0就表示该浮点数是正数.

其次再来讨论M的存储，前面提到了规约浮点数的1<=M<2, 也就是说M可以表示为1.x1x2x3…. 这样的形式.其中 .x1x2x3… 是小数部分. 对于规约浮点数，这个1是一定会存在的. 所以IEEE754规定，在计算机内保存规约浮点数的M时， 默认这个数的第一位为 1，因此可以被舍去，这样子就可以节省一位有效数字位，使得32（64）位规约浮点数可以保存24（53）位的有效数字。

最后再来讨论E的存储. E的存储和解释稍微麻烦一点. 首先，E是一个无符号整数（unsign int ）,着意味着当E为8位时，其取值范围为0到255；若E为11位其取值范围为0到2047。但是我们知道，科学计数法中的E可以是负数，因此，E的真实值必须减去一个中间值。对于8位的E应减去127，对于11位的E应减去1023；所以对于单精度浮点数而言, 它的指数最多到达255-127=128（其实后面我们会知道255其实是保留值，表示±inf或者NaN），所以单精度浮点的范围大概是 (-2^128, 2^128)，即（-3.4e38, 3.4e38）, 而双精度浮点数的范围大概是(-2^1024, 2^1024), 即（-1.79e308 , 1.79e308）.

比如说，2^9保存为32位浮点数时, 其E必须保存为9+127=136=10001000

而还原E的真实值分成以下三种情况

  1. E不全为0也不全为1，这时浮点数就是规约浮点数. 则E直接减去127就是原先规约浮点数在十进制下的指数值.
  2. E全为0。这时，浮点数就是非规约浮点数. 这时浮点数的指数E为1-127=-126（双精度下是1-1023=-1022），有效数字M不再加上第一位1，而是依旧保留为0.x1x2x3… ，如果此时M也全是0的话，则就表示±0（正还是负取决于符号位S）
  3. E全为1。这时如果有效数字M全为0，则表示±inf（到底是+inf还是-inf取决于符号位s）；如果有效数字M不全为0，表示这个数是一个NaN

等等, 我们前面只说了将规约浮点数转换为二进制的时候E就是+127，没讲非规约浮点数转换为二进制该怎么转换的呢~ 但其实非规约浮点数我们一般用不到. 因为非规约浮点数太小了. 你想啊，E都是0了. 则真实值就是-127了.

下面来完成176.0625作为单精度浮点存储的二进制表示，上面以及将它变成二进制，即

(176.0625)= (10110000.0001)

再将10110000.0001转换为IEEE754的数学模型（注意M=1.01100000001在[1,2)范围内）

10110000.0001 = 1.01100000001x2^7

将M的1去掉，得到01100000001,再将01100000001扩展为23位尾数，即01100000001000000000000

而E是7+127=134=10000110（无符号整数）

最后，将符号位（S）、阶码（E）与尾数（M）合并起来，最终得到176.0625的表示形式

0 10000110 01100000001000000000000

下面说说为什么需要非规约浮点数. 即E为0的浮点数. 这是因为规约浮点数有空隙!

这怎么讲? 规约浮点数绝对值的最小值是1.0*2^(-126), 即M=1，E=1(因为规约浮点数不允许E全为0或者1， 1-127=-126)，即首先规约浮点数无法表示真正的0.

而浮点数绝对值的次小值是 (1.0+2^(-23))*2^(-126), 其次，绝对值最小的正浮点数A=1.0*2^(-126)到0的距离是A到比他稍微大一点的相邻浮点数B= (1.0+2^(-23))*2^(-126)的距离的2^23倍!!! 一图胜千言

即 |B-A| = |A|*2^23. 这就很不合理~ 出现断崖式的递减到0（而不是均匀的递减到0）. 不过通过上图我们也知道，实数在计算机中的存储也不是连续的啦~

而出现了非规约浮点的话，就解决了这个问题。首先0的表示上面已经说了，0是一个非规约浮点数. 其次非规约浮点数填补了[0,A]之间的空白。非规约浮点数的指数是1-127=-126. 即 2^(-126), 而最小的正的非规约浮点数是

2^(-126)*2^(-23)=2^(-149), 最大的非规约浮点数是 2^(-126)^(1-2^(-23))=2^(-126)-2^(-149)

所以非规约浮点之间的间隙是 2^(-149), 而且均匀的填满了 [0,A], 甚至最大的非规约和最小的规约之间的距离也是2^(-149), 即

这就完美的解决了仅有规约浮点导致的断崖式递减为0的问题

好了，说了浮点数在计算机中的存储结构, 下面看看一道C语言的题目

1
2
3
4
5
6
7
8
9
10
11
12
13
14
#include <stdio.h>
int main()
{
    int num = 9;
    float *pFloat = &num;
    printf("num 的值为：%d\n",num);
    printf("*pFloat 的值为：%f\n",*pFloat);

    *pFloat = 9.0;
    printf("num 的值为：%d\n",num);
    printf("*pFloat 的值为：%f\n",*pFloat);
    
    return 0;
}

运行结果是

这是为什么呢?

首先第一行是没有问题的. 因为num本身就是一个整型. 而将9赋给一个单精度浮点之后为什么变成0了呢? 这里首先要明确一点——数据在计算机中宏统统是二进制，计算机并不知道这些个字节表示的到底是浮点还是整型还是其他什么类型的数据. num在计算机中的存储为

00000000000000000000000000001001

而再把这个二进制串赋给一个float之后. 计算机就会依照float的规则来读取并使用这个二进制（即计算机中的二进制到底是什么，取决于用途，而非存储形式）.

而按照非规约浮点表示法，S=0，E=00000000（全部是0），所以E还原回去就是1-127=-126 最后

00000000000000000000000000001001的M是 00000000000000000001001

综上，这个32位二进制作为float被读取其实表示

(-1)^0*0.00000000000000000001001*2^(-126)=1.001*2^(-146)（一个非常小的正数）

当然，上面的表达式是二进制的.

而转成十进制的话，因为float精度最多到6~7位, 所以是0

而9.0用int表示却是十进制的1092567616 ，这又是为什么呢? 还是那句话——计算机中的二进制到底是什么，取决于用途，而非存储形式

所以我们首先将9.0在计算机中的存储二进制搞出来.

首先：浮点数9.0的二进制表示为1001.0, 即1.001*2^3, 符号位S=0

其次, 有效数字M=100 0000 0000 0000 0000 0000（共23位（100后加上20个0）其中最高位1默认被省略）

最后，指数E=3+127=130，即E=10000010

综上，9.0作为单精度浮点在计算机中的32位存储二进制为

01000001000100000000000000000000

而你现在命令计算机以int的格式去读取它，不就是十进制的1091567616吗?

最后我们来谈谈浮点数精度丢失的问题. 首先你得知道double和float的精度——其实就是M的有效数字个数

float：2^23 = 8388608（其实应该是2^(-23), 但是从2^23也可以估计出2^-23的小数点位数），一共七位，这意味着最多能有7位有效数字，但绝对能保证的为6位，也即float的精度为6~7位有效数字；

double：2^52 = 4503599627370496，一共16位，同理，double的精度为15~16位有效数字。

也就是说8位数的int和float之间转换就会出精度丢失的问题

下面进一步说明这个问题

Google的首席科学家Vincent Vanhoucke在Udacity上开设的Deep Learning课程上讲述了如下一个例子, 计算

理论值显然应该为1。

但是下面的程序给出了

1
2
3
4
5
6
7
8
9
10
int main() {
	double eps1 = 1e9, eps2=1e-6;
	for (int i = 0;i<1e6; i++)
	{
		eps1+=eps2;
	}
	eps1-=1e9;
	printf("%lf\n", eps1);
	return 0;
}

这样的答案~ 为什么会有这个误差呢? 这就不得不用到我们刚刚说过的浮点数精度问题.

一个超大的数和一个超小的数相加，会发生什么呢?

假设当前的浮点数变量可以保存6位有效数字的数值。那么，数值12345与1.234相加的理论值应该是12346.234。但由于当前只能保存6位有效数值，这个正确的理论值会被截断为12346.2，这就出现了0.034的误差。当有很多这样的大数与小数相加时，截断误差就会逐步累积，导致最后的计算结果出现大的偏差。

不要小看这些误差，Google工程师的这个例子可以说明，在100万次求和后，截断误差的累积量已经非常可观了。

一般的，在计算机求和的过程中，一个大数和小数的相加会因为浮点数的有限精度，而导致截断误差的出现。所以，当需要对一系列的数值做加法时，一个好的技巧是将这些数由大到小做排列，再逐个相加，而避免两个相差极大的数直接相加. 而如果应用场景一定要做出这样的大数与小数的求和，怎么解决精度丢失问题呢?

一个直观的想法就是——你这个误差不就是大数吃小数导致的嘛~ 我只要在相加的时候将大数和小数分别拆成A和B两部分, 确保2个A之间相加不会丢失精度，也确保2个B部分之间相加不会丢失精度，随着计算次数不断的增多, B部分的累加也很可观了以至于再和A相加也不会丢失精度了. 则再将B的部分融入A的部分. 这个想法就是Kahan求和， Kahan求和的目的就是降低浮点精度丢失.

用此思想，我们不难改进上面的程序

1
2
3
4
5
6
7
8
9
10
11
12
13
int main() {
	double sum = 1e9, c = 0, eps = 1e-6,s,t;
	for (int i = 0; i<1e6;i++)
	{
		s = eps+c; // 每次要加的数(1e-6)补上之前因为精度丢失而导致的量,这里不会丢失精度
		t = sum+s; // sum每次加上的是s而不再是eps, 这里将导致精度丢失, 因为是一个超大的数和一个超小的数相加, 可以理解为t比真实值小
		c = s-(t-sum); // t比sum增加的量比真实应该增加的量s少,少的量恰好就是丢失的精度(用c保存当前次丢失的精度然后在下一次补回去就是kahan算法的精髓), 这里的运算不会丢失精度,因为不是超大的数和超小的数
		sum = t;
	}
	sum-=1e9;
	printf("%lf\n", sum);
	return 0;
}

我们欣喜的发现

即使用Kahan求和能够得到近乎于理论的精度数值.

关键就在于将补偿量c记录在了不会丢失精度的代码中.

参考

【1】https://blog.csdn.net/Marvine/article/details/89817146

【2】http://blog.sina.com.cn/s/blog_5fb3f1250100xodv.html

【3】https://blog.csdn.net/jvandc/article/details/81176294

【4】https://blog.csdn.net/zhangpinghao/article/details/8138732

【5】https://www.jianshu.com/p/988020aa2a5b