Gibt es ein Problem mit Gleitkommaberechnungen?

Question

Bedenken Sie den folgenden Code:

0,1 + 0,2 == 0,3 ->

0,1 + 0,2 ->

Warum treten diese Ungenauigkeiten auf?

Answer 1

硬件设计师的视角

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

1. 概述

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

2. 标准

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

在IEEE-754标准中，硬件设计师可以允许任何误差/ε值，只要它小于最后一位的半个单位，并且结果仅对于一个操作必须小于最后一位的半个单位。这解释了为什么在重复操作时，误差会累积。对于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
• C99 hexfloat notation：0x1.999999999999ap-4

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

• 十进制：0.1
• 类似于C99 hexfloat notation的表示：0x1.99999999999999...p-4，其中...表示无尽的9的序列。

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

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

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

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

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

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

在实践中，这个精度问题意味着你需要使用舍入函数将浮点数舍入到你感兴趣的小数位数，然后再显示它们。

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

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

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

其中abs是绝对值函数。myToleranceValue需要根据你的特定应用选择，它与你准备允许多少“摇摆余地”以及你将要比较的最大数字有很大关系（由于精度损失问题）。要注意你所选择的语言中的“epsilon”风格常量。这些常量可以用作容差值，但它们的有效性取决于你正在处理的数字的大小，因为与大数计算可能超过epsilon阈值。