浮点数运算是否存在问题？

Question

考虑以下代码：

0.1 + 0.2 == 0.3  ->  false

0.1 + 0.2         ->  0.30000000000000004

为什么会出现这些不准确性？

Answer 1

硬件设计师的视角

我认为我应该从硬件设计师的角度来添加一些观点，因为我设计和构建浮点硬件。了解错误的起源可能有助于理解软件中发生的情况，最终，我希望这有助于解释浮点错误发生的原因，并且似乎随着时间的推移而累积。

1. 概述

从工程角度来看，由于执行浮点计算的硬件只需要在最后一位的一半以下具有一个单位的误差，所以大多数浮点运算都会有一定的误差。因此，大多数硬件只会停止在只需要产生一个单操作中最后一位的一半以下误差的精度，这在浮点除法中尤为棘手。什么构成一个单操作取决于单元接受的操作数数量。对于大多数单元来说，是两个，但有些单元接受3个或更多的操作数。因此，无法保证重复的操作会产生理想的误差，因为这些误差会随着时间的推移而累积。

2. 标准

大多数处理器遵循IEEE-754标准，但有些使用非规范化或不同的标准。例如，IEEE-754中有一种非规范化模式，允许以牺牲精度的方式表示非常小的浮点数。然而，下面将介绍IEEE-754的规范化模式，这是典型的操作模式。

在IEEE-754标准中，硬件设计师可以选择任何误差/epsilon值，只要它小于最后一位的一半以下单位，并且结果只需要在一个操作中小于最后一位的一半以下单位。这解释了为什么在重复的操作中，误差会累积。对于IEEE-754双精度来说，这是第54位，因为有53位用于表示浮点数的数值部分（规范化部分），也称为尾数（例如5.3e5中的5.3）。接下来的几节将更详细地介绍各种浮点运算中硬件错误的原因。

3. 除法舍入误差的原因

浮点除法中误差的主要原因是用于计算商的除法算法。大多数计算机系统使用乘法的逆来计算除法，主要是在Z=X/Y和Z = X * (1/Y)中。除法是迭代计算的，即每个周期计算一些商的位数，直到达到所需的精度，对于IEEE-754来说，这是任何误差小于最后一位的一半以下的内容。Y的倒数表（1/Y）称为慢除法中的商选择表（QST），商选择表的位数通常是基数的宽度，或者每次迭代计算的商的位数加上几个保护位。对于IEEE-754标准的双精度（64位），它将是除法器的基数大小加上几个保护位k，其中k>=2。因此，例如，一个每次计算2位商（基数为4）的典型商选择表将是2+2= 4位（加上几个可选位）。

3.1 除法舍入误差：倒数的近似

商选择表中的倒数取决于除法方法：慢除法（如SRT除法）或快速除法（如Goldschmidt除法）；每个条目根据除法算法进行修改，以尽量减小误差。无论如何，所有倒数都是实际倒数的近似，并引入一定的误差。慢除法和快速除法方法都是迭代计算商，即每个步骤计算一定数量的商位数，然后从被除数中减去结果，除数重复这些步骤，直到误差小于最后一位的一半以下。慢除法方法在每个步骤中计算一定数量的商位数，通常更便宜，而快速除法方法在每个步骤中计算可变数量的位数，通常更昂贵。除法方法最重要的部分是它们大多数都依赖于通过近似的倒数进行重复乘法，因此它们容易出错。

4. 其他操作中的舍入误差：截断

所有操作中舍入误差的另一个原因是IEEE-754允许的不同截断模式。有截断、向零舍入、向最近舍入（默认）、向下舍入和向上舍入。所有方法都会为单个操作引入小于最后一位的一半单位的误差。随着时间和重复操作的进行，截断也会累积到结果误差中。这种截断误差在涉及某种形式的重复乘法的指数运算中尤为棘手。

5. 重复操作

由于执行浮点计算的硬件只需要在单个操作中产生一个单位的最后一位的一半以下误差的结果，如果不加以观察，随着重复操作，误差将会增加。这就是为什么在需要有界误差的计算中，数学家使用方法，例如使用IEEE-754的最近舍入的偶数位，因为随着时间的推移，错误更有可能相互抵消，以及区间算术结合变种的IEEE 754舍入模式来预测舍入误差并进行修正。由于相对误差较低，与其他舍入模式相比，最近舍入的偶数位（在最后一位）是IEEE-754的默认舍入模式。

请注意，默认舍入模式，最近舍入的偶数位，保证一个操作的误差小于最后一位的一半单位。仅使用截断、向上舍入和向下舍入可能导致大于最后一位的一半单位但小于最后一位单位的误差，因此除非在区间算

Answer 2

二进制浮点数运算就是这样的。在大多数编程语言中，它基于IEEE 754标准。问题的关键在于，数字在这种格式中表示为一个整数乘以二的幂；分母不是二的幂的有理数（例如0.1，即1/10）无法被精确表示。

对于标准的binary64格式中的0.1，其表示可以精确地写为

• 0.1000000000000000055511151231257827021181583404541015625（十进制），或
• 0x1.999999999999ap-4（C99十六进制浮点表示法）。

相比之下，有理数0.1，即1/10，可以精确地写为

• 0.1（十进制），或
• 0x1.99999999999999...p-4（类似于C99十六进制浮点表示法，其中...表示无尽的9的序列）。

您程序中的常量0.2和0.3也将是它们真实值的近似值。恰好0.2最接近的double大于有理数0.2，但最接近的double小于有理数0.3。0.1和0.2的和最终大于有理数0.3，因此与代码中的常量不一致。

关于浮点数运算问题的一个相当全面的处理是《计算机科学家应该了解的浮点运算知识》。更易理解的解释，请参见floating-point-gui.de。

普通的十进制（基数10）数字也存在相同的问题，这就是为什么像1/3这样的数字最终会变成0.333333333...的原因。

你刚好遇到了一个（3/10）在十进制系统中很容易表示，但在二进制系统中无法表示的数字。这种情况也是相反的（在某种程度上）：1/16在十进制中是一个丑陋的数字（0.0625），但在二进制中看起来像十进制中的1/10000那样整洁（0.0001）** - 如果我们习惯于在日常生活中使用二进制数字系统，你甚至会看到那个数字并本能地理解你可以通过不断折半来得到它。

当然，这并不完全是浮点数在内存中存储的方式（它们使用一种科学计数法的形式）。然而，它确实说明了一个问题，即二进制浮点精度错误往往会出现，因为我们通常感兴趣的“现实世界”数字往往是十的幂次数 - 但只是因为我们日常使用的是十进制数系统。这也是为什么我们会说71%而不是“每7个中的5个”（71%是一个近似值，因为任何十进制数都无法精确地表示5/7）。

所以不，二进制浮点数并没有损坏，它们只是像其他任何基数N的数系统一样不完美 :)

在实践中，这种精度问题意味着您需要使用舍入函数将浮点数舍入到您感兴趣的小数位数之前再显示。

您还需要用允许一定容差的比较来替换等式测试，这意味着：

不要使用if (x == y) { ... }

而是使用if (abs(x - y) < myToleranceValue) { ... }。

其中abs是绝对值函数。需要根据您的特定应用选择myToleranceValue - 这与您准备允许多少“摆动空间”以及您要比较的最大数字有很大关系（由于精度损失问题）。请注意您所选择的语言中的“epsilon”样式常量。这些常量可以用作容差值，但其有效性取决于您正在使用的数字的大小（由于大数计算可能超过epsilon阈值）。