8

我只能假设这是一个错误。第一个断言通过，而第二个断言失败：

double sum_1 =  4.0 + 6.3;
assert(sum_1 == 4.0 + 6.3);

double t1 = 4.0, t2 = 6.3;

double sum_2 =  t1 + t2;
assert(sum_2 == t1 + t2);

如果不是错误，为什么？

4

7 回答 7

13

您正在比较浮点数。不要那样做，浮点数在某些情况下具有固有的精度误差。相反，取两个值之差的绝对值，并断言该值小于某个小数 (epsilon)。

void CompareFloats( double d1, double d2, double epsilon )
{
    assert( abs( d1 - d2 ) < epsilon );
}

这与编译器无关，与浮点数的实现方式有关。这是 IEEE 规范：

http://www.eecs.berkeley.edu/~wkahan/ieee754status/IEEE754.PDF

于 2009-08-26T23:06:36.343 回答
12

这也是让我感到痛苦的事情。

是的，由于舍入误差，永远不应该比较浮点数是否相等，您可能知道这一点。

但在这种情况下，您正在计算t1+t2，然后再次计算它。当然，这必须产生相同的结果？

这是可能发生的事情。我敢打赌你是在 x86 CPU 上运行的，对吗？x86 FPU 使用 80 位作为其内部寄存器，但内存中的值存储为 64 位双精度数。

所以t1+t2首先以 80 位精度计算，然后——我想——以 64 位精度存储到内存中sum_2——然后进行一些舍入。对于断言，它被加载回浮点寄存器，并t1+t2再次以 80 位精度计算。因此，现在您正在比较sum_2以前四舍五入为 64 位浮点值的 与t1+t2以更高精度（80 位）计算的 ——这就是为什么这些值不完全相同的原因。

编辑那么为什么第一个测试通过了？在这种情况下，编译器可能会在编译时进行评估4.0+6.3并将其存储为 64 位的数量——既用于赋值，也用于断言。因此，正在比较相同的值，并且断言通过。

第二次编辑这是为代码的第二部分（gcc，x86）生成的汇编代码，带有注释——几乎遵循上面概述的场景：

// t1 = 4.0
fldl    LC3
fstpl   -16(%ebp)

// t2 = 6.3
fldl    LC4
fstpl   -24(%ebp)

// sum_2 =  t1+t2
fldl    -16(%ebp)
faddl   -24(%ebp)
fstpl   -32(%ebp)

// Compute t1+t2 again
fldl    -16(%ebp)
faddl   -24(%ebp)

// Load sum_2 from memory and compare
fldl    -32(%ebp)
fxch    %st(1)
fucompp

有趣的旁注：这是在没有优化的情况下编译的。使用 编译时-O3，编译器会优化所有代码。

于 2009-08-26T23:16:37.350 回答
3

我已经在我的 Intel Core 2 Duo 上复制了您的问题，并查看了汇编代码。这是发生了什么：当你的编译器评估时t1 + t2，它确实

load t1 into an 80-bit register
load t2 into an 80-bit register
compute the 80-bit sum

当它存储到sum_2它时

round the 80-bit sum to a 64-bit number and store it

然后==比较将 80 位和与 64 位和进行比较，它们是不同的，主要是因为小数部分 0.3 不能使用二进制浮点数精确表示，所以您正在比较“重复小数”（实际上重复二进制）已被截断为两种不同的长度。

真正令人恼火的是，如果您使用gcc -O1or进行编译gcc -O2，gcc 在编译时会执行错误的算术运算，问题就会消失。根据标准，这可能是可以的，但这只是 gcc 不是我最喜欢的编译器的另一个原因。

PS 当我说==将 80 位和与 64 位和进行比较时，当然我的意思是它比较了 64 位和的扩展版本。你最好想想

sum_2 == t1 + t2

解决为

extend(sum_2) == extend(t1) + extend(t2)

和

sum_2 = t1 + t2

解决为

sum_2 = round(extend(t1) + extend(t2))

欢迎来到浮点的奇妙世界！

于 2009-08-26T23:46:37.587 回答
3

在比较浮点数的接近度时，您通常希望测量它们的相对差异，其定义为

if (abs(x) != 0 || abs(y) != 0) 
   rel_diff (x, y) = abs((x - y) / max(abs(x),abs(y))
else
   rel_diff(x,y) = max(abs(x),abs(y))

例如，

rel_diff(1.12345, 1.12367)   = 0.000195787019
rel_diff(112345.0, 112367.0) = 0.000195787019
rel_diff(112345E100, 112367E100) = 0.000195787019

这个想法是测量数字共有的前导有效数字的数量；如果你取 0.000195787019 的 -log10，你会得到 3.70821611，这大约是所有示例共有的前 10 位数字的数量。

如果您需要确定两个浮点数是否相等，您应该执行类似的操作

if (rel_diff(x,y) < error_factor * machine_epsilon()) then
  print "equal\n";

其中 machine epsilon 是可以在所使用的浮点硬件的尾数中保存的最小数字。大多数计算机语言都有一个函数调用来获取这个值。error_factor 应该基于您认为在计算数字 x 和 y 时舍入错误（和其他）会消耗的有效位数。例如，如果我知道 x 和 y 是大约 1000 次求和的结果，并且不知道求和数字的任何界限，我会将 error_factor 设置为大约 100。

试图将这些添加为链接，但由于这是我的第一篇文章，所以不能：

  • en.wikipedia.org/wiki/Relative_difference
  • en.wikipedia.org/wiki/Machine_epsilon
  • en.wikipedia.org/wiki/Significand（尾数）
  • en.wikipedia.org/wiki/Rounding_error
于 2009-08-27T02:09:46.840 回答
2

可能在其中一种情况下，您最终会将 64 位双精度寄存器与 80 位内部寄存器进行比较。查看 GCC 针对这两种情况发出的汇编指令可能会很有启发性……

于 2009-08-26T23:12:28.700 回答
1

双精度数的比较本质上是不准确的。例如，您经常会发现0.0 == 0.0返回false。这是由于 FPU 存储和跟踪数字的方式。

维基百科说：

测试平等是有问题的。两个数学上相等的计算序列很可能产生不同的浮点值。

您将需要使用 delta 来为您的比较提供一个容差，而不是一个确切的值。

于 2009-08-26T23:08:13.467 回答
0

这个“问题”可以通过使用这些选项来“修复”：

-msse2 -mfpmath=sse

如本页所述：

http://www.network-theory.co.uk/docs/gccintro/gccintro_70.html

一旦我使用了这些选项，两个断言都通过了。

于 2009-08-27T21:52:38.620 回答