使用双打时,为什么 (x / (y * z)) 与 (x / y / z) 不同?

When using doubles, why isn't (x / (y * z)) the same as (x / y / z)?

这部分是学术性的,就我的目的而言,我只需要四舍五入到小数点后两位;但我很想知道是什么导致了两个略有不同的结果。

这是我为将其缩小到最简单的实现而编写的测试:

@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

任何人都可以详细解释导致 1.000000000000000X 处的值不同的原因吗?

我在答案中寻找的一些要点是: 精度在哪里丢失? 首选哪种方法,为什么? 哪个才是正确的? (在纯数学中,两者不可能都是对的。也许两者都错了?) 这些算术运算有没有更好的解决方案或方法?

当然,操作顺序与 double 不精确 :

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

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

那是因为双除法经常会导致精度损失。所述损失可能因除法顺序而异。

除以 7d 时,实际结果已经失去了一些精度。然后只有你将错误的结果除以 60.

除以7d * 60时,只需使用一次除法,因此只损失一次精度。

请注意,双倍乘法有时也会失败,但这不太常见。

它与 double 类型的实现方式以及浮点类型不提供与其他更简单的数字类型相同的精度保证这一事实有关。尽管以下答案更具体地与求和有关,但它也通过解释浮点数学运算中如何不能保证无限精度来回答您的问题:Why does changing the sum order returns a different result?. Essentially you should never attempt to determine the equality of floating-point values without specifying an acceptable margin of error. Google's Guava library includes DoubleMath.fuzzyEquals(double, double, double) to determine the equality of two double values within a certain precision. If you wish to read up on the specifics of floating-point equality this site is quite useful; the same site also explains floating-point rounding errors。总而言之:您的计算的预期值和实际值不同,因为计算之间的舍入因运算顺序而异。

我看到一堆问题告诉您如何解决这个问题,但除了 "floating-point roundoff error is bad, m'kay?" 之外没有一个真正解释发生了什么,所以让我试一试。首先让我指出,此答案中没有任何内容特定于 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 所做的那样。每当您执行的算术运算将精度位数增加到超出数据类型可以表示的范围时,它就会再次 发生。这些效果中的任何一个都称为 roundoff error。有几种不同的方式可以发生这种情况:

  • 加法和减法:如果您要添加的值中的一个与另一个具有不同的指数,您将得到额外的精度位数,如果它们足够多,则最不重要的将需要被删除。例如,2.718 和 121.0 都是可以在 float6d 数据类型中精确表示的值。但是,如果您尝试将它们加在一起:

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

    四舍五入为 1.237 x 10^2,即 123.7,精度下降两位数。

  • 乘法:结果的位数大约为两个操作数的位数之和。这将产生一些舍入误差,如果您的操作数已经有很多有效数字。例如,121 x 2.718 给你

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

    四舍五入为 3.289 x 10^2 或 328.9,再次降低两位数的精度。

    但是,请记住,如果您的操作数是 "nice" 数字,没有很多有效数字,浮点格式可能可以准确表示结果,因此您不必处理舍入误差。例如,2.3 x 140 给出

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

    没有舍入问题。

  • 分裂:这就是事情变得混乱的地方。除法几乎 总是 会导致一些舍入误差,除非你除以的数字恰好是基数的幂(在这种情况下,除法只是一个数字移位,或二进制位移)。例如,取两个非常简单的数字 3 和 7,将它们相除,得到

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

    可以表示为 float6d 的最接近该数字的值是 4.286 x 10^-1 或 0.4286,这与确切结果明显不同。

正如我们将在下一节中看到的那样,舍入引入的误差会随着您执行的每项操作而增加。所以 如果你正在使用 "nice" 数字,就像你的例子一样,通常最好尽可能晚地进行除法运算 因为这些是最有可能引入的操作舍入错误进入您的程序,其中 none 之前存在。

舍入误差分析

一般来说,如果你不能假设你的数字是"nice",舍入误差可以是正的也可以是负的,并且很难仅仅根据操作来预测它会走向哪个方向。这取决于所涉及的具体值。查看 2.718 z 的舍入误差图作为 z 的函数(仍然使用 float6d 数据类型):

实际上,当您处理使用数据类型的完全精度的值时,通常更容易将舍入误差视为随机误差。查看该图,您可能会猜到误差的大小取决于运算结果的数量级。在这种特殊情况下,当 z 的顺序为 10-1 时,2.718 z 的顺序也为 10-1,因此它将是 0.XXXX 形式的数字。最大舍入误差是最后一位精度的一半;在这种情况下,"the last digit of precision" 是指 0.0001,因此舍入误差在 -0.00005 和 +0.00005 之间变化。在2.718 z跳升到下一个数量级的点,即1/2.718 = 0.3679,可以看到舍入误差也跳升了一个数量级。

您可以使用众所周知的 techniques of error analysis 来分析某个量级的随机(或不可预测的)错误如何影响您的结果。具体来说,对于乘法或除法,结果中的 "average" 相对误差可以通过在每个操作数 求积 中添加相对误差来近似计算 - 即对它们进行平方,将它们相加,然后取平方根。对于我们的 float6d 数据类型,相对误差在 0.0005(对于 0.101 这样的值)和 0.00005(对于 0.995 这样的值)之间变化。

让我们将 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 * yz / x中的平均相对误差为0.0001414,z / x * y中的相对误差为

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

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

血淋淋的细节

您可能对您在问题中提出的具体计算感到好奇,而不仅仅是平均值。对于该分析,让我们切换到二进制算术的真实世界。大多数系统和语言中的浮点数使用 IEEE standard 754. For 64-bit numbers, the format 表示,指定 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

您可以使用 this converter 或 Internet 上的任何其他工具来使用这些值。

当您计算第一个表达式时,450/(7*60),处理器首先进行乘法运算,得到 420,或者

1.1010010000000000000000000000000000000000000000000000 x 2^00000001000

然后用 450 除以 420。得到 15/14,即

1.0001001001001001001001001001001001001001001001001001001001001001001001...

二进制。现在,the Java language specification 表示

Inexact results must be rounded to the representable value nearest to the infinitely precise result; if the two nearest representable values are equally near, the one with its least significant bit zero is chosen. This is the IEEE 754 standard's default rounding mode known as round to nearest.

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

1.0001001001001001001001001001001001001001001001001001 x 2^00000000000

这大约是十进制的 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)恰好是准确值的最准确表示。这是一个例子,说明离开分裂直到最后对你有好处。但同样,只要您使用的值不使用数据类型的完整精度,这条规则就成立。一旦开始使用不精确(四舍五入)的值,您就不再通过先进行乘法来保护自己免受进一步的舍入错误。

让我们稍微简化一下。你想知道的是为什么 450d / 420450d / 7 / 60(特别是)给出不同的结果。

让我们看看除法在IEE双精度浮点格式中是如何进行的。无需深入实施细节,它基本上是 XOR-ing 符号位,从被除数的指数中减去除数的指数,除以尾数,并对结果进行归一化。

首先,我们应该以 double 的正确格式表示我们的数字:

450    is  0 10000000111 1100001000000000000000000000000000000000000000000000

420    is  0 10000000111 1010010000000000000000000000000000000000000000000000

7      is  0 10000000001 1100000000000000000000000000000000000000000000000000

60     is  0 10000000100 1110000000000000000000000000000000000000000000000000

我们先450除以420

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

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

看起来不错,现在是尾数:

1.1100001000000000000000000000000000000000000000000000 / 1.1010010000000000000000000000000000000000000000000000 == 1.1100001 / 1.101001。有几种不同的方法可以做到这一点,我稍后会详细讨论它们。结果是1.0(001)(可以验证here)。

现在我们应该规范化结果。让我们看看 guard、round 和 sticky 位值:

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

保护位已设置,圆形和粘性位未设置。我们四舍五入到最接近的值(IEEE 的默认模式),我们被困在我们可以四舍五入的两个可能值之间。由于 lsb 是 0,我们添加 1。这给了我们四舍五入的尾数:

0000000100100100100100100100100100100100100100100101

结果是

0 10000000101 0000000100100100100100100100100100100100100100100101

十进制表示为 64.28571428571429

现在我们必须将它除以 60...但是您已经知道我们失去了一些精度。 450 除以 420 根本不需要四舍五入,但在这里,我们已经不得不将结果 至少四舍五入一次 。但是,为了完整起见,让我们完成这项工作:

64.28571428571429 除以 60

符号位=0

指数 = 10000000101b - 10000000100b + 01111111111b == 01111111110b

尾数 = 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

Tl;博士:

第一师给了我们:

0 01111111111 0001001001001001001001001001001001001001001001001001

最后一个部门给了我们:

0 01111111111 0001001001001001001001001001001001001001001001001010

区别仅在于最后 2 位,但我们可能会损失更多 - 毕竟,要获得第二个结果,我们必须舍入 两次而不是 none!

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

IEEE 长除法规定的方式(here 是一些很好的例子;它基本上是常规的长除法,但用二进制而不是十进制),而且速度很慢。这就是你的电脑所做的。

还有一个更快但不太准确的选项,乘以逆。先求除数的倒数,再做乘法