诗海

我们越谦卑，就离真理越近

0%

详解浮点数

发表于 更新于 分类于
本文字数： 7.8k 阅读时长 ≈ 7 分钟

这篇文章跟我的上一篇"使用牛顿法解决开平方问题"直接相关，当时在书写牛顿法的迭代代码时，需要判断循环终止的条件。Princeton大学的版本使用的是while (Math.abs(t - a / t) > EPSILON * t)这样的语句，初看貌不惊人，实际上思考起来却大有玩味。不要看就这简简单单的几行程序，我相信会有很多人在第一次写的时候，会像我一样写成while (Math.abs(a - t * t) > EPSILON)或者是比较当前值和下一次运算迭代值的while (Math.abs(a / t - t) / 2.0 > EPSILON)这个样子。实际上，后两种的写法都是有问题的。那么这里面的学问是什么呢？这就需要我们对于计算机的浮点数运算有深入的理解才能进行回答。

我们知道，数学中的小数是无穷无尽的，但是计算机只能用有限的位数（单精度浮点数32位，双精度浮点数64位）进行表达。这就意味着，一定有一些数是计算机无法准确表达的。然而不幸的是，我们无论是在机器学习，计算机图形学，数据分析等等领域，用到浮点数的计算又相当之广泛，这就给计算机的浮点数表示形式带来了很大的挑战。如果不清楚浮点数的本质原理，就会很容易在程序编写中造成意想不到的问题。战战兢兢，如履薄冰，这是我们在使用浮点数时应有的谨慎态度。

浮点数的表达形式和精度

浮点数的表达形式遵循IEEE 754中的规范，即使用以2为基数的有限小数位的科学记数法（Scientific Notation）来表示和存储。以32位浮点数为例，表示形式为：

  • 1位：符号位（sign），0表示正数，1表示负数
  • 8位：指数位（exponent），范围0-255
  • 23位：小数位（fraction），通常表示1.xxx小数点后面的有效数字，即\(0.f_1f_2f_3...f_{23}\)，也有例外

对于浮点数的数值表示，可以用以下这个公式进行描述： \[ V = (-1)^S \times M \times 2^E \]

根据exponent的不同，我们将浮点数分为三种类型：

  1. 规格化的值

    当Exponent不全为0且不全为1时，这时是浮点数的一般情况，此时\(M = 1 + fraction = 1.f_1f_2f_3...f_{23}\)，\(E = exponent - Bias\)，这里\(Bias = 2^{8 - 1} - 1 = 127\)。规格化值在\([1 \times 2^{-126}, (2-2^{-23}) \times 2^{127}]\)区间分布，越远离数轴0值分布越稀疏，越接近数轴0值分布越密集（只考虑非负数，负浮点数与正浮点数以原点0在数轴上对称分布）。

  2. 非规格化的值

    当Exponent全为0时，此时属于非规格化的浮点数，\(M = fraction = 0.f_1f_2f_3...f_{23}\)，\(E = 1 - Bias = -126\)，采用这样的定义方式是为了均匀的表示非常接近于0的数，保证在\([0, 2^{-127})\)区间内也有浮点数的分布（全按规格化值来的话这一区间是没有浮点数数值的），所有的非规格化值在\([0, (2 - 2^{-23}) \times 2^{-126}]\)区间上均匀分布（只考虑非负数）。

  3. 特殊值

    当Exponent全为1时，此时分为两种情况：若fraction部分不全为0，则表示特殊值“NaN”（Not a Number）；若fraction部分全为0，则表示正无穷（\(+\infty\)）或负无穷（\(-\infty\)）。

Type Sign Actual Exp Exp Exponent field Fraction field Value
零 0 -126 0 0000 0000 000 0000 0000 0000 0000 0000 0.0
负零 1 -126 0 0000 0000 000 0000 0000 0000 0000 0000 -0.0
一 0 0 127 0111 1111 000 0000 0000 0000 0000 0000 1.0
负一 1 0 127 0111 1111 000 0000 0000 0000 0000 0000 -1.0
最小非规格化数 * -126 0 0000 0000 000 0000 0000 0000 0000 0001 \(\pm 2^{-149} \approx \pm 1.4 \times 10^{-45}\)
中间非规格化数 * -126 0 0000 0000 100 0000 0000 0000 0000 0000 \(\pm 2^{-127} \approx \pm 5.88 \times 10^{-39}\)
最大非规格化数 * -126 0 0000 0000 111 1111 1111 1111 1111 1111 \(\pm (1 - 2^{-23}) \times 2^{-126} \approx \pm 1.18 \times 10^{-38}\)
最小规格化数 * -126 1 0000 0001 000 0000 0000 0000 0000 0000 \(\pm 2^{-126} \approx \pm 1.18 \times 10^{-38}\)
最大规格化数 * 127 254 1111 1110 111 1111 1111 1111 1111 1111 \(\pm (2 - 2^{-23}) \times 2^{127} \approx \pm 3.4 \times 10^{38}\)
正无穷 0 128 255 1111 1111 000 0000 0000 0000 0000 0000 \(+ \infty\)
负无穷 1 128 255 1111 1111 000 0000 0000 0000 0000 0000 \(- \infty\)
不是一个数 * 128 255 1111 1111 non zero NaN

64位的浮点数和32位的差不多，表示形式为：

  • 1位：符号位（sign）
  • 11位：指数位（exponent）
  • 52位：小数位（fraction）

类比32位的浮点数，可以知道\(Bias = 2^{11 - 1} - 1 = 1023\)，64位浮点数的表示范围是：

非规格化浮点数\([0, (2 - 2^{-52}) \times 2^{-1022}]\)(十进制表示约为\([0.0, 2.2 \times 10^{-308}]\))，

规格化浮点数\([1 \times 2^{-1022}, (2 - 2^{-52}) \times 2^{1023}]\)(十进制表示约为\([2.2 \times 10^{-308}, 1.8 \times 10^{308}]\))

以下的图表描述了单精度浮点数（32位）和双精度浮点数（64位）二进制和对应的近似十进制转化：

浮点数的运算和舍入

不管是使用32位还是64位，都是使用有限的位数来表达数学中的小数，这就意味着一定有一些小数是无法准确使用计算机浮点数进行表达的。最典型的例子例如：使用浮点数表示0.7

根据计算可以得知，十进制的0.7转化为二进制的表示是：0.1011(0011)* = \(1.011(0011) \times 2^{-1}\)，因此转化为32位浮点数是Sign = 0, Exp = 126, frac = 011(0011)*，32位浮点数写为: 0011 1111 0011 0011 0011 0011 0011 0011，使用十六进制简写表示为：0x3f333333. 此时这个32位浮点数表示的实际值，实际上是: 0.699999988079071044921875

除了0.7之外，还有一些数是无法用浮点数准确表示的，比如以0.1为代表的所有以10为基数的负指数幂（1e-1, 1e-15, 3e-5这样的十进制科学计数。这个背后的问题是因为二进制无法在有限位准确表示1/5，就像我们说十进制中无法在有限位准确表示1/3一样），这也就是为什么在涉及要求精确的金融数据计算时，不能使用IEEE 754浮点数的缘故，一般来说，推荐使用BigDecimal, int, long等类型进行金融计算。

以上在将十进制的0.7转化为二进制数时，已经使用到了浮点数的舍入（也就是保留二进制小数点右边前23位）：浮点数的舍入规则分为向偶数舍入、向零舍入、向上舍入、向下舍入四种情况，一般默认采用向偶数舍入(round-to-even, 也被成为向最接近的值舍入round-to-nearest)。例如：

当我们希望保留二进制小数点右边前2位时，我们会将\(10.00011_2(2\frac{3}{32})\)向下舍入为\(10.00_2(2)\)，\(10.0110_2(2\frac{3}{16})\)向下舍入为\(10.01_2(2\frac{1}{4})\)，这些值不是正中间值，因此直接进行舍入。我们将\(10.11100_2(2\frac{7}{8})\)向上舍入为\(11.00_2(3)\)，而将\(10.10100_2(2\frac{5}{8})\)向下舍入为\(10.10_2(2\frac{1}{2})\)，这些值都是两个可能值的中间值，此时根据向偶数舍入原则（使最低有效位为零）来操作。

在浮点数进行运算的过程中，同样涉及到了舍入，而且可能会产生意想不到的效果。例如：

使用32位浮点数计算表达式\((3.14 + 1e10) - 1e10\)求值得到0.0，因为舍入，这里的3.14会丢失；相反，计算表达式\(3.14 + (1e10 - 1e10)\)则会得到值3.14。这里，

1e10用二进制表示为：0101 0000 0001 0101 0000 0010 1111 1001(0x501502f9) = 10'000'000'000

比1e10小的上一个数：0101 0000 0001 0101 0000 0010 1111 1000(0x501502f8) = 9'999'998'976

比1e10大的下一个数：0101 0000 0001 0101 0000 0010 1111 1010(0x501502fa) = 10'000'001'024

可以看到，在32位浮点数中，1e10和上一个浮点数，以及下一个浮点数的距离之差高达1024！因此我们也就可以理解，当1e10先与3.14进行加法运算后，所得到的float值肯定是向下舍入为1e10，然后再与1e10相减变为0.0，在这一过程中，3.14就被舍入掉了。

（相比较而言，在64位浮点数表示下，1e10为：0x4202A05F20000000 = 10'000'000'000， 它的上一个数是：0x4202A05F1FFFFFFF = 9'999'999'999.999'998'092'651'367'187'5，下一个数是：0x4202A05F20000001 = 10'000'000'000.000'001'907'348'632'812'5，距离之差已经非常小了！）

由此，我们也就可以明白，32位浮点数采用23位fraction和64位浮点数采用52位fraction的精度差距在那里了：

对于32位float，可以准确地表示\([-2^{24}, 2^{24}]\)之间（相当于十进制数\([-16777216, 16777216]\)）的任何整数，对于再往数轴两侧延伸的整数，则只能挑有限的进行表示（比如1e10可以准确表示，但是1e10 + 1就不能表示）

对于64位double，可以准确地表示\([-2^{53}, 2^{53}]\)之间的任何整数，对于再往数轴两侧延伸的整数，则只能挑有限的进行表示。

非规格化浮点数

前面以及说过，当exponent部分取最小值0时，则目前的浮点数就是非规格化的（Denormal number）。非规格化数表示那些非常接近0的数，它是为了更好的处理算术下溢（underflow）而引入的。算术下溢是说浮点数计算的结果小于当前浮点数类型能表示的正常浮点数的最小值（最小非规格化数），此时即发生下溢。在没有非规格化浮点数的老式设计中，发生下溢的数会被直接刷新成0，而使用非规格化浮点数则可以使用一个最接近真实值的非规格化浮点数来表示运算结果，减少一定的精度损失。

非规格化浮点数在靠近0的地方均匀分布，这点与规格化浮点数越向数轴两侧延伸则两个相邻值距离越大的特点不同。

一般来说，不推荐在程序中使用非规格化的浮点数进行加、减、乘、除的运算。这是因为计算机硬件的缘故，会造成对非规格化浮点数的计算不能由处理器直接进行，而是需要进行陷入和特殊指令处理。一般来说，计算机对于非规格化浮点数的计算速度会比规格化浮点数慢10-100倍左右。

浮点数的比较

在谈及浮点数的比较之前，我们首先引入一个概念，叫做ULP（unit in the last place）。从浮点数分布来说，它描述的是任意一个给定的浮点数与比它fraction大1的下一个浮点数的距离值。正如上文所说，在32位float中，1e10的ULP就是1024.0。在Java中，可以使用Math.ulp()查看当前浮点数的ULP值。

浮点数的分布特点，可以这样来表述：在两个相邻的2的指数幂之间（比如\([2.0, 4.0]\)之间，或是\([2^{-126}, 2^{-125}]\)之间）的所有浮点数，它们的ULP相同；越向数轴两侧延伸，则ULP以2的倍数增长。唯一的例外是0和非规格化数，这些数相邻之间的距离都是相等的，以32位浮点数为例，ULP = \(2 ^ {-149}\)(也就是最小非规格化的浮点数值)

有1个非常特殊的ULP值，那就是1.0时的ULP。在C语言的<float.h>中，这一特殊值被定义为FLT_EPSILON和DBL_EPSILON，即最小的浮点数值使得1.0 + x != 1.0：

在32位float中，这个值是\(2^{-23}\)，也即是1.1920928955078125E-7

在64位double中，这个值是\(2^{-52}\)，也即是2.22044604925031308084726333618E-16

为什么说这个值很特殊呢？这是因为在接下来要介绍的浮点数比较方法中，采用相对距离比较时EPSILON的取值一般默认就是这个值。另外，这个值与最小规格化浮点数的乘积刚好是最小非规格化浮点数，这意味着在两个规格化浮点数的比较过程中，判断条件右侧的表达式不会出现舍入为0.0的情况。

比较两个浮点数a, b是否相近主要分为3种方法：

  • 绝对值比较（absolute difference），也就是我们说的abs(a - b)。这种方法问题很大，正如上文讨论ULP时提到的，使用Math.abs(a - b) <= EPSILON进行比较时，若a，b的浮点数很大，则ULP很大，此时一个很小的EPSILON值无法准确的衡量两个数是否接近，甚至只有a，b相等的情况才会满足不等式条件。这也就解释了为什么引言中提到的两种绝对值比较的方法，都是错误的。除非我们可以保证if (Math.abs(a - b) <= EPSILON)中的EPSILON能够大于1ULP，允许相近值的存在，否则原式就等价于if (Math.abs(a - b) == 0.0)

  • 相对距离比较（relative distance），这个是目前主流的比较两个数是否接近的方法，它比较的是一个比例，当两个浮点数之间的差值所占比例很小时，就认为两个浮点数是相近的。（这里分母取两浮点数中较大的、较小的、中间值都是可以的，以下仅以较大值为例）Math.abs(a - b) / Math.max(a, b) <= EPSILON一般写作if (Math.abs(a - b) <= Math.max(a, b) * EPSILON)我个人感觉这是因为由于浮点数本质是科学记数法，所以乘法的运算会比除法更快一些？这里EPSILON取值就有学问了，当取1.0时的ULP特殊值时，不等式的成立等同于要么两浮点数相等，要么相邻（距离不超过1ULP）(特殊情况：在2的指数幂和它左侧的浮点数比较时，容忍范围扩大为原先的2倍)。当然我们也可以将EPSILON取值为1.0时ULP特殊值的若干倍，这样就是可以在两个浮点数相邻若干个ULP时，都认为它们是接近的，相当于扩大了"相近"的容忍范围。

  • 编辑距离比较（edit distance），也就是说在a, b两个浮点数之间还有多少个有效的浮点数。使用这种方法的函数比较难以计算，一般只用于已知两个距离非常接近的数之间，还有就是两个非规格化浮点数的比较情况。想象一个极端的情况: 32位浮点数的最小非规格化数0x00000001（\(2^{-149}\)）和它相邻的下一个浮点数0x0000002（\(2^{-148}\)），两者之间只相差了一个ULP，但是差值占比却是50%！在这样的情况下，理论上使用相对距离比较就是不行的；在实际操作上，由于max * epsilon的值小于最小非规格化浮点数，会直接舍入为0.0。因此，两个非规格化的浮点数的近似比较一般不使用相对距离。

对引言的回答

在有了前面所有的基础之后，我们再来看引言中提到的循环的控制条件while (Math.abs(t - a / t) > EPSILON * t)

首先，由于t最后收敛于\(\sqrt a\)，因此t和a / t的double值必然会越来越接近。在这两个值中，t是较大的那一个（前文证明过数列的单调性），因此循环条件的反面就是if (Math.abs(t - a / t) <= EPSILON * Math.max(t, a / t))，这就是标准的相对距离比较时用到的条件判断式。

其次，我们再来说EPSILON，这里EPSILON定义为1e-15。浮点数是不能准确表示1e-15的，这里实际上是它的近似值1.00000000000000007770539987666E-15(0x3CD203AF9EE75616)。想必已经发现了，它是一个大于1.0时ULP的一个浮点数。这就意味着，我们允许两个浮点数相距4.5个ULP。如果使用1e-16或者比1.0时ULP更小的浮点数，那么条件判断式就退化为必须两个浮点数相等了，此时就存在一种无限循环下去的风险。

最后，我们说这个循环条件在面对所有规格化的double值时表现都是没有问题的，它的缺陷在于如果尝试进行两个非规格化的浮点数的距离比较，就会出现max * epsilon下溢为0的情况。但是，这里的好处是哪怕是最小的非规格化浮点数Double.MIN_VALUE（\(2^{-1074}\)），它的Sqrt值也是一个规格化的浮点数(\(2^{-537} > 2^{-1022}\))，此外，在将t初始化为a的条件下，第一次的t赋值语句t = (t + a / t) / 2.0实际上就是(a + 1.0) / 2.0，这里a直接被舍入掉了，就如同（float条件下）1e10 + 3.14一样。因此，之后的t值一直是在规格化浮点数的范围内进行改变的。

参考资料