当使用双精度浮点数时,为什么(x / (y * z))不等同于(x / y / z)?

24

这在一定程度上是学术性的,因为就我的目的而言,我只需要将其四舍五入到小数点后两位;但我很想知道是什么导致了两个稍微不同的结果。

这是我编写的测试,以将其缩小到最简实现:

@Test
public void shouldEqual() {
  double expected = 450.00d / (7d * 60);  // 1.0714285714285714
  double actual = 450.00d / 7d / 60;      // 1.0714285714285716

  assertThat(actual).isEqualTo(expected);
}

但它的输出结果是失败的:

org.junit.ComparisonFailure: 
Expected :1.0714285714285714
Actual   :1.0714285714285716

有人可以详细解释一下,在底层发生了什么导致X为1.000000000000000时的值不同吗?

我希望回答中包括以下几点: 精度损失在哪里? 哪种方法更受推荐,为什么? 哪一个是真正正确的?(在纯数学中,两者都不可能是正确的。也许两者都错了?) 这些算术操作是否有更好的解决方案或方法?


1
@UwePlonus 我不这么认为。那个问题及其答案是关于如何消除影响的,而不是对底层情况的真正解释。 - David Z
2
因为它实际上会进行不同的计算。第一行给出了 450.00d / (420d)(在计算 7d * 60 的第一步中没有精度损失)。第二行首先计算 450.00d / 7d 并存储结果,由于计算机存储浮点数的方式,这一步会有少量精度损失,然后将该结果除以 60。您可以在此处了解浮点数的工作原理:http://floating-point-gui.de/formats/fp/ - Nils O
1
@UwePlonus 不,我认为这并不符合Ben在这里所要求的详细解释。 - David Z
3
如果你对Java如何进行浮点数运算的奇特现象感兴趣,那么你可能需要阅读我的文章存档,链接在这里http://blogs.msdn.com/b/ericlippert/archive/tags/floating+point+arithmetic/ 和 http://ericlippert.com/tag/floating-point-arithmetic/。虽然这些文章是关于C#和JavaScript的,但其中大部分内容同样适用于Java。 - Eric Lippert
1
即使是小数,也存在精度限制的问题。例如,使用5个有效数字进行计算,8/3/2 = 2.6667/2 = 1.3334,而8/(3*2) = 8/6 = 1.3333。 - phuclv
显示剩余7条评论
5个回答

43
我看到很多问题告诉你如何解决这个问题,但没有一个真正解释了发生了什么,除了“浮点舍入误差很糟糕”。所以让我来试一下。首先,让我指出,这个答案中没有任何特定于Java的内容。舍入误差是固定精度表示数字所固有的问题,因此在C等语言中也会出现相同的问题。

十进制数据类型中的舍入误差

作为一个简化的例子,想象一下我们有一种计算机,它本地使用无符号的十进制数据类型,让我们称其为float6d。数据类型的长度为6位数字:4位用于尾数,2位用于指数。例如,数字3.142可以表示为

3.142 x 10^0

这将被存储为6位数字

503142

前两位是指数加50,后四位是尾数。这种数据类型可以表示从0.001 x 10^-509.999 x 10^+49之间的任何数字。

实际上,那并不是真的。如果你想要表示3.141592?或者3.1412034?或者3.141488906?很遗憾,这种数据类型不能存储超过四位精度的数字,所以编译器必须将具有更多数字的任何内容舍入以适应数据类型的限制。如果你写

float6d x = 3.141592;
float6d y = 3.1412034;
float6d z = 3.141488906;

然后编译器将这三个值转换为相同的内部表示形式,即3.142 x 10^0(记住,它存储为503142),因此x == y == z成立。

关键在于,有一整个实数范围映射到相同的底层数字序列(或在实际计算机中为位)。具体而言,任何满足3.1415 <= x <= 3.1425(假设使用“半偶数”舍入)的x都会被转换为表示为503142以在内存中存储。

每次程序将浮点数值存储在内存中时,都会进行此舍入。第一次发生的时间是当您在源代码中写入常量时,就像我上面使用xyz一样。每当执行增加数据类型所能表示的精度位数的算术运算时,它都会再次发生。这两种效果之一称为舍入误差。这可能会以几种不同的方式发生:

  • Addition and subtraction: if one of the values you're adding has a different exponent from the other, you will wind up with extra digits of precision, and if there are enough of them, the least significant ones will need to be dropped. For example, 2.718 and 121.0 are both values that can be exactly represented in the float6d data type. But if you try to add them together:

       1.210     x 10^2
    +  0.02718   x 10^2
    -------------------
       1.23718   x 10^2
    

    which gets rounded off to 1.237 x 10^2, or 123.7, dropping two digits of precision.

  • Multiplication: the number of digits in the result is approximately the sum of the number of digits in the two operands. This will produce some amount of roundoff error, if your operands already have many significant digits. For example, 121 x 2.718 gives you

       1.210     x 10^2
    x  0.02718   x 10^2
    -------------------
       3.28878   x 10^2
    

    which gets rounded off to 3.289 x 10^2, or 328.9, again dropping two digits of precision.

    However, it's useful to keep in mind that, if your operands are "nice" numbers, without many significant digits, the floating-point format can probably represent the result exactly, so you don't have to deal with roundoff error. For example, 2.3 x 140 gives

       1.40      x 10^2
    x  0.23      x 10^2
    -------------------
       3.22      x 10^2
    

    which has no roundoff problems.

  • Division: this is where things get messy. Division will pretty much always result in some amount of roundoff error unless the number you're dividing by happens to be a power of the base (in which case the division is just a digit shift, or bit shift in binary). As an example, take two very simple numbers, 3 and 7, divide them, and you get

       3.                x 10^0
    /  7.                x 10^0
    ----------------------------
       0.428571428571... x 10^0
    

    The closest value to this number which can be represented as a float6d is 4.286 x 10^-1, or 0.4286, which distinctly differs from the exact result.

作为下一节中我们将看到的,由四舍五入引入的误差随着每次操作而增加。因此,如果您使用的是“好”数字(如您的示例),通常最好尽可能晚地进行除法运算,因为这些运算最有可能在程序中引入舍入误差。
舍入误差分析
一般来说,如果您不能假设您的数字是“好”的,则舍入误差可能是正的或负的,并且很难根据操作预测它将朝哪个方向发展。它取决于涉及的具体值。看一下这个以float6d数据类型为例的2.718 z的舍入误差随着z变化的图表:

roundoff error for multiplication by 2.718

实际上,当您使用数据类型的全精度值时,将舍入误差视为随机误差通常更容易。从图中可以看出,您可能能够猜测误差的大小取决于操作结果的数量级。在这种特定情况下,当 z 的数量级为 10-1 时,2.718 z 也是数量级为 10-1 的数字,因此它将成为形式为 0.XXXX 的数字。最大舍入误差然后是精度的最后一位数字的一半;在这种情况下,我指的是 0.0001 的“精度的最后一位数字”,因此舍入误差在 -0.00005 和 +0.00005 之间变化。在 2.718 z 跳到下一个数量级(即 1/2.718 = 0.3679)的点上,您可以看到舍入误差也会跳一个数量级。
你可以使用众所周知的 误差分析技术 来分析某个大小的随机(或不可预测)误差如何影响你的结果。具体而言,对于乘法或除法,你的结果中的“平均”相对误差可以通过将每个操作数的相对误差 平方 - 即将它们平方、相加并取平方根来近似计算。使用我们的 float6d 数据类型,相对误差在 0.0005(例如 0.101)和 0.00005(例如 0.995)之间变化。

relative error in values between 0.1 and 1

让我们将0.0001作为值xy的相对误差的粗略平均值。则x * yx / y的相对误差由以下公式给出:

sqrt(0.0001^2 + 0.0001^2) = 0.0001414

这个因子比每个单独值的相对误差大sqrt(2)倍。

当涉及到多个操作时,您可以为每个浮点运算应用此公式多次。例如,对于z / (x * y)x * y的相对误差平均为0.0001414(在此十进制示例中),然后z / (x * y)的相对误差为

sqrt(0.0001^2 + 0.0001414^2) = 0.0001732

注意,平均相对误差随每次操作增加而增加,具体地说,它随你执行的乘除运算数量的平方根增长。
同样地,在 z / x * y 中,z / x 的平均相对误差为 0.0001414,z / x * y 的相对误差为。
sqrt(0.0001414^2 + 0.0001^2) = 0.0001732

所以,在这种情况下是一样的。这意味着对于任意值,平均而言,两个表达式引入的误差大约相同。(理论上是这样。我见过这些操作在实践中表现非常不同,但那是另一回事。)

详细信息

您可能对问题中提出的具体计算感到好奇,而不仅仅是一个平均值。为了进行分析,让我们转到二进制算术的真实世界。在大多数系统和语言中,浮点数使用IEEE标准754进行表示。对于64位数字,格式指定52位专用于尾数,11位专用于指数,1位专用于符号。换句话说,当以2为底写成浮点数时,其形式为:

1.1100000000000000000000000000000000000000000000000000 x 2^00000000010
                       52 bits                             11 bits

前导的1没有明确存储,它构成了第53位。此外,您应该注意,用于表示指数的11位实际上是真实指数加1023。例如,这个特定的值是7,即1.75 x 22。尾数是二进制下的1.75,或者1.11,指数是1023 + 2 = 1025,二进制下为10000000001,因此存储在内存中的内容为

01000000000111100000000000000000000000000000000000000000000000000
 ^          ^
 exponent   mantissa

但这并不重要。 你的例子还涉及450。
1.1100001000000000000000000000000000000000000000000000 x 2^00000001000

和60,

1.1110000000000000000000000000000000000000000000000000 x 2^00000000101

你可以使用此转换器或互联网上的其他许多转换器来尝试这些值。
当您计算第一个表达式450/(7*60)时,处理器首先进行乘法运算,得到420。
1.1010010000000000000000000000000000000000000000000000 x 2^00000001000

然后它将450除以420,得到15/14。
1.0001001001001001001001001001001001001001001001001001001001001001001001...

在二进制下。现在,Java语言规范指出:

不精确的结果必须舍入为最接近无穷精确结果的可表示值;如果两个最接近的可表示值同样接近,则选择其最低有效位为零的值。这是IEEE 754标准的默认舍入模式,称为四舍五入至最近。

在64位IEEE 754格式中,15/14的最近可表示值为

1.0001001001001001001001001001001001001001001001001001 x 2^00000000000

这段文字涉及编程相关内容。第一段文字中的数字 1.0714285714285714 在十进制下约等于 1.0714285714285714,更准确地说,这是最不精确的十进制值,可以唯一确定该特定二进制表示法。另一方面,如果首先计算 450 / 7,则结果为 64.2857142857...,或用二进制表示为:
1000000.01001001001001001001001001001001001001001001001001001001001001001...

最接近的可表示值为

1.0000000100100100100100100100100100100100100100100101 x 2^00000000110

这是64.28571428571429180465... 注意到二进制尾数的最后一位(与精确值相比)由于舍入误差而发生了变化。将其除以60即可得到

1.000100100100100100100100100100100100100100100100100110011001100110011...

看一下结尾:模式不同了!重复的是0011,而不是另一种情况下的001。最接近可表示的值为

1.0001001001001001001001001001001001001001001001001010 x 2^00000000000

这与最后两位的其他运算顺序不同:它们是10而不是01。十进制等效值为1.0714285714285716。

如果您查看确切的二进制值,应该可以清楚地看到导致此差异的特定舍入:

1.0001001001001001001001001001001001001001001001001001001001001001001001...
1.0001001001001001001001001001001001001001001001001001100110011001100110...
                                                     ^ last bit of mantissa

在这种情况下,前面的结果数值上为15/14恰好是最准确的表示。这是一个将除法留到最后的好处的例子。但是,只要你处理的值没有使用数据类型的完整精度,这个规则才能保持有效。一旦你开始使用不精确(四舍五入)的值,通过先进行乘法运算就不能再保护自己免受进一步的舍入误差的影响了。

5
这与double类型的实现方式有关,以及浮点类型不能像其他更简单的数字类型一样提供相同的精度保证。虽然下面的答案更具体地讨论了求和问题,但它也通过解释浮点数数学运算没有无限精度的保证来回答了你的问题:为什么改变求和顺序会返回不同的结果?。实际上,在没有指定可接受的误差范围的情况下,你永远不应该尝试确定浮点值的相等性。Google的Guava库包括DoubleMath.fuzzyEquals(double, double, double),用于在一定精度内确定两个double值的相等性。如果你想了解浮点相等性的具体细节,这个网站非常有用;同样的网站还解释了浮点舍入误差。总之,你的计算期望值和实际值之间的差异是由于操作顺序导致的舍入不同。

4

让我们简化一下。你想知道的是为什么 450d / 420450d / 7 / 60(具体来说)会得到不同的结果。

让我们看看在IEEE双精度浮点格式中如何执行除法。不深入实现细节,基本上是将除数的指数从被除数的指数中减去、除以尾数并对结果进行归一化。

首先,我们应该用适当的格式表示我们的数字 double:

450    is  0 10000000111 1100001000000000000000000000000000000000000000000000

420    is  0 10000000111 1010010000000000000000000000000000000000000000000000

7      is  0 10000000001 1100000000000000000000000000000000000000000000000000

60     is  0 10000000100 1110000000000000000000000000000000000000000000000000

首先让我们将 450 除以 420

首先是符号位,它是 0 (0 xor 0 == 0)。

然后是指数。 10000000111b - 10000000111b + 1023 == 10000000111b - 10000000111b + 01111111111b == 01111111111b

看起来很好,现在是尾数:

1.1100001000000000000000000000000000000000000000000000 / 1.1010010000000000000000000000000000000000000000000000 == 1.1100001 / 1.101001。有几种不同的方法可以做到这一点,稍后我会谈一些关于它们的事情。结果是 1.0(001) (您可以在 这里 验证)。

现在我们应该标准化结果。让我们看一下保护、舍入和粘性位值:

0001001001001001001001001001001001001001001001001001 0 0 1

保护位为0,我们不进行任何舍入。结果以二进制表示为:

0 01111111111 0001001001001001001001001001001001001001001001001001

在十进制中表示为1.0714285714285714

现在让我们类比地将450除以7

符号位=0

指数=10000000111b - 10000000001b + 01111111111b == -01111111001b + 01111111111b + 01111111111b == 10000000101b

尾数=1.1100001 / 1.11 == 1.00000(001)

舍入:

0000000100100100100100100100100100100100100100100100 1 0 0

Guard位被设置,舍入和粘滞位未被设置。我们采用 IEEE 的默认舍入方式——四舍五入,并且正好处于两个可能的取值之间。由于最低有效位是 0,我们加上 1。这给了我们舍入后的尾数: 0000000100100100100100100100100100100100100100100101 结果为 0 10000000101 0000000100100100100100100100100100100100100100100101 在十进制中表示为 64.28571428571429
现在我们需要将其除以 60……但你已经知道我们失去了一些精度。将 450 除以 420 根本不需要舍入,但在这里,我们已经至少舍入了一次结果。但是,为了完整起见,让我们完成这项工作: 64.28571428571429 除以 60 符号位 = 0 Exponent = 10000000101b - 10000000100b + 01111111111b == 01111111110b

Mantissa = 1.0000000100100100100100100100100100100100100100100101 / 1.111 == 0.10001001001001001001001001001001001001001001001001001100110011

圆整和移位:

0.1000100100100100100100100100100100100100100100100100 1 1 0 0

1.0001001001001001001001001001001001001001001001001001 1 0 0

舍入与之前的情况一样,我们得到尾数:0001001001001001001001001001001001001001001001001010
当我们将其向左移动1位时,我们将其加到指数中,得到
指数 = 01111111111b 因此,结果为: 0 01111111111 0001001001001001001001001001001001001001001001001010 在十进制中表示为1.0714285714285716
简而言之:
第一次除法给出了: 0 01111111111 0001001001001001001001001001001001001001001001001001 最后一次除法给出了: 0 01111111111 0001001001001001001001001001001001001001001001001010 唯一的区别在于最后的2位,但我们可能会失去更多-毕竟,要获得第二个结果,我们必须舍入两次而不是零次!

现在,关于尾数除法。浮点除法有两种主要实现方式。

IEEE长除法规定的方法(这里有一些很好的例子;它基本上是常规的长除法,但使用的是二进制而不是十进制),速度相当慢。这就是您的计算机所做的。

还有一种更快但精度较低的选择,即倒数乘法。首先找到除数的倒数,然后进行乘法。


1
由于双重除法经常导致精度丢失。该损失取决于除法的顺序。
当你用 7d 进行除法时,实际结果已经失去了一些精度。然后,您只需将错误的结果除以 60
当您通过 7d * 60 进行除法时,您只需要使用一次除法,因此仅会失去一次精度。
请注意,双重乘法有时也会出现问题,但这种情况要少得多。

请注意,双倍乘法有时也会失败,但这种情况要少得多,只适用于整数参数。对于非整数,这种情况同样普遍;例如,0.1*0.1 != 0.01 - user2357112

0

当然,操作的顺序与 双精度浮点数不精确 这一事实混合在一起:

450.00d / (7d * 60) --> a = 7d * 60 --> result = 450.00d / a

对比

450.00d / 7d / 60 --> a = 450.00d /7d --> result = a / 60

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接