浮点比较不可重复性

Question

我和我的博士。学生在物理数据分析环境中遇到了一个问题，我可以对此进行一些见解。我们有代码可以分析来自大型强子对撞机实验之一的数据，该实验给出了不可重复的结果。特别是，从在同一台机器上运行的相同二进制文件获得的计算结果在连续执行之间可能会有所不同。我们知道不可再现性的许多不同来源，但排除了常见的嫌疑人。

我们将问题归结为在比较名义上具有相同值的两个数字时（双精度）浮点比较运算的不可再现性。由于分析中的先前步骤，这种情况有时可能会发生。我们刚刚找到一个测试数字是否小于 0.3 的示例（请注意，我们从不测试浮点值之间的相等性）。事实证明，由于探测器的几何形状，计算有时可能会产生恰好为 0.3（或其最接近的双精度表示形式）的结果。

我们非常清楚比较浮点数的陷阱，以及 FPU 中精度过高可能影响比较结果的可能性。我想回答的问题是“为什么结果不可重现？”是否因为 FPU 寄存器加载或其他 FPU 指令没有清除多余的位，因此先前计算的“剩余”位影响了结果？ （这似乎不太可能）我在另一个论坛上看到一个建议，进程或线程之间的上下文切换也可能导致浮点比较结果发生变化，因为 FPU 的内容存储在堆栈上，因此被截断。对这些=或其他可能的解释的任何评论将不胜感激。

Answer 1

据猜测，发生的情况是，您的计算通常在 FPU 内部以一些额外的精度进行，并且仅在特定点进行舍入（例如，当您将结果分配给值时）。

然而，当发生上下文切换时，必须保存和恢复 FPU 的状态 - 并且至少很有可能在上下文切换中不保存和恢复这些额外的位。当这种情况发生时，这可能不会导致重大变化，但如果（例如）您稍后从每个金额中减去固定金额并乘以剩下的金额，则差异也会成倍增加。

需要明确的是：我怀疑“剩余”的部分会是罪魁祸首。相反，它会丢失额外的位，从而导致计算中略有不同的点进行舍入。

Answer 2

什么平台？

大多数 FPU 可以在内部存储比 IEEE 双精度表示更高的精度 - 以避免中间结果中的舍入误差。通常编译器会切换到速度/准确性 - 请参阅 http ://msdn.microsoft.com/en-us/library/e7s85ffb(VS.80).aspx

Answer 3

程序是多线程的吗？

如果是，我会怀疑竞争条件。

如果不是，程序执行是确定性的。在给定相同输入的情况下获得不同结果的最可能原因是未定义的行为，即程序中的错误。读取未初始化的变量、过时的指针、覆盖堆栈上某些 FP 编号的最低位等。可能性是无限的。如果您在 Linux 上运行它，请尝试在 valgrind 下运行它，看看它是否发现了一些错误。

顺便说一句，你是如何将问题缩小到 FP 比较的？

（远景：硬件故障？例如，RAM 芯片故障可能会导致数据在不同情况下读取方式不同。不过，这可能会很快使操作系统崩溃。）

任何其他解释都是令人难以置信的——操作系统或硬件中的错误不会很长时间不被发现。

Answer 4

我做了这个：

#include <stdio.h>
#include <stdlib.h>

typedef long double ldbl;

ldbl x[1<<20];

void hexdump( void* p, int N ) {
  for( int i=0; i<N; i++ ) printf( "%02X", ((unsigned char*)p)[i] );
}

int main( int argc, char** argv ) {

  printf( "sizeof(long double)=%i\n", sizeof(ldbl) );

  if( argc<2 ) return 1;

  int i;
  ldbl a = ldbl(1)/atoi(argv[1]);

  for( i=0; i<sizeof(x)/sizeof(x[0]); i++ ) x[i]=a;

  while(1) {
    for( i=0; i<sizeof(x)/sizeof(x[0]); i++ ) if( x[i]!=a ) {
      hexdump( &a, sizeof(a) );
      printf( " " );
      hexdump( &x[i], sizeof(x[i]) );
      printf( "\n" );
    }
  }

}

使用 /Qlong_double 使用 IntelC 进行编译，因此它产生了这个：

;;;     for( i=0; i<sizeof(x)/sizeof(x[0]); i++ ) if( x[i]!=a ) {

        xor       ebx, ebx                                      ;25.10
                                ; LOE ebx f1
.B1.9:                          ; Preds .B1.19 .B1.8
        mov       esi, ebx                                      ;25.47
        shl       esi, 4                                        ;25.47
        fld       TBYTE PTR [?x@@3PA_TA+esi]                    ;25.51
        fucomp                                                  ;25.57
        fnstsw    ax                                            ;25.57
        sahf                                                    ;25.57
        jp        .B1.10        ; Prob 0%                       ;25.57
        je        .B1.19        ; Prob 79%                      ;25.57
[...]
.B1.19:                         ; Preds .B1.18 .B1.9
        inc       ebx                                           ;25.41
        cmp       ebx, 1048576                                  ;25.17
        jb        .B1.9         ; Prob 82%                      ;25.17

并启动了 10 个具有不同“种子”的实例。正如你所看到的，它比较了
10 字节长从内存中加倍，其中一个位于 FPU 堆栈上，因此在这种情况下
操作系统不能保持完整的精度，我们肯定会看到错误。
好吧，他们仍然在运行而没有检测到任何东西......这并不是真的
令人惊讶的是，因为 x86 有立即保存/恢复整个 FPU 状态的命令，
无论如何，不​​能保持完全精度的操作系统将被完全破坏。

所以要么是一些独特的操作系统/CPU/编译器，要么是不同的比较结果
在更改程序中的某些内容并重新编译它或其或其
程序中的错误，例如。缓冲区溢出。

Answer 5

CPU 的内部 FPU 可以以比 double 或 float 更高的精度存储浮点。每当寄存器中的值必须存储在其他地方时，包括当内存被换出到缓存中时（我知道这一事实）以及该核心上的上下文切换或操作系统中断听起来像是另一个简单的源，这些值都必须进行转换。当然，操作系统中断或上下文切换或非热内存交换的时间是应用程序完全无法预测和控制的。

当然，这取决于平台，但您的描述听起来像是在现代桌面或服务器（因此 x86）上运行。

Answer 6

我将合并 David Rodriguez 和 Bo Persson 的一些评论并进行大胆的猜测。

使用SSE3指令时会不会发生任务切换？基于此关于使用 SSE3 指令的英特尔文章 保存寄存器状态的命令 FSAVE 和 FRESTOR 已被 FXSAVE 和 FXRESTOR 取代，它们应该处理
累加器的全长。

在 x64 机器上，我认为“不正确”的指令可能包含在某些外部编译库中。

Answer 7

您肯定遇到了 GCC Bug n°323，与其他点一样out 是由于 FPU 精度过高造成的。

解决方案：

• 执行计算
• 使用 SSE（或 AVX，这是 2016 年...）使用 -ffloat-store 编译开关 。来自海湾合作委员会文档。

不要将浮点变量存储在寄存器中，并禁止其他可能改变浮点值是从寄存器还是内存中获取的选项。
此选项可防止在 68000 等机器上出现不需要的过高精度，其中（68881 的）浮动寄存器保持比双精度更高的精度。
x86 架构也是如此。
对于大多数程序来说，过高的精度只会带来好处，但也有少数程序依赖于 IEEE 浮点的精确定义。
在修改此类程序以将所有相关的中间计算存储到变量中之后，请对此类程序使用 -ffloat-store 。