我和我的博士学生在物理数据分析环境中遇到了一个问题,我可以使用一些洞察力。我们的代码可以分析来自其中一个LHC实验的数据,这些实验可以产生不可重复的结果。特别是从中得到的计算结果 相同的二进制,跑上去 同一台机器 连续执行可能会有所不同。我们知道许多不同的不可再生性来源,但已排除了通常的嫌疑人。
当比较两个名义上具有相同值的数字时,我们已经将问题跟踪到(双精度)浮点比较操作的不可再现性。由于分析中的先前步骤,偶尔会发生这种情况。我们刚刚找到一个例子来测试一个数字是否小于0.3(注意我们从不测试浮点值之间的相等性)。事实证明,由于探测器的几何形状,计算有时可能产生精确到0.3(或其最接近的双精度表示)的结果。
我们非常清楚比较浮点数以及FPU中过度精度可能影响比较结果的缺陷。我想回答的问题是“为什么结果不可复制?”是因为FPU寄存器加载或其他FPU指令没有清除多余的位,因此先前计算中的“剩余”位会影响结果吗? (这似乎不太可能)我在另一个论坛上看到一个建议,即进程或线程之间的上下文切换也可能导致浮点比较结果的变化,因为FPU的内容存储在堆栈中,因此被截断。任何关于这些=或其他可能的解释的评论将不胜感激。
猜测,正在发生的事情是,您的计算通常在FPU内部执行一些额外的精度,并且仅在特定点处舍入(例如,当您将结果赋值给某个值时)。
然而,当存在上下文切换时,必须保存和恢复FPU的状态 - 并且这些额外位是至少非常公平的。 不 在上下文切换中保存和恢复。当它发生时,这可能不会引起重大变化,但是如果(例如)你稍后从每个变量中减去固定数量并乘以剩下的数量,那么差异也会成倍增加。
需要明确的是:我怀疑“遗留下来”是罪魁祸首。相反,它会丢失额外的位,导致计算中稍微不同点的舍入。
什么平台?
大多数FPU可以在内部存储比ieee双重表示更高的精度 - 以避免中间结果中的舍入误差。通常有一个编译器切换来交换速度/准确性 - 请参阅 http://msdn.microsoft.com/en-us/library/e7s85ffb(VS.80).aspx
我做的:
#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状态的命令,
无论如何,一个不能保持完全精度的操作系统将完全被破坏。
所以要么是它的一些独特的OS / cpu /编译器,要么是不同的比较结果
在更改程序中的某些内容并重新编译它之后生成它
程序中的错误,例如。缓冲区溢出。
该程序是多线程的吗?
如果是的话,我会怀疑是一场竞争。
如果不是,程序执行是确定性的。在给定相同输入的情况下获得不同结果的最可能的原因是未定义的行为,即程序中的错误。读取未初始化的变量,过时的指针,覆盖堆栈上某些FP编号的最低位等。可能性是无穷无尽的。如果你在linux上运行它,尝试在valgrind下运行它,看看它是否揭示了一些错误。
顺便说一下,你是如何将问题缩小到FP比较的?
(远射:硬件故障?例如,RAM芯片出现故障可能会导致数据在不同情况下以不同方式被读取。但是,这可能会很快崩溃操作系统。)
任何其他解释都是不可信的 - 操作系统或硬件中的错误不会被发现很长时间。
CPU的内部FPU可以以高于double或float的精度存储浮点。只要寄存器中的值必须存储在其他地方,包括当存储器被换出到缓存中时,这些值就必须被转换(我知道这是事实),并且该核心上的上下文切换或OS中断听起来像另一个简单的来源。当然,OS中断或上下文切换的时序或非热存储器的交换是完全不可预测的,并且应用程序无法控制。
当然,这取决于平台,但您的描述听起来像是在现代桌面或服务器上运行(因此x86)。
我将合并David Rodriguez和Bo Persson的一些评论并做出疯狂猜测。
在使用SSE3指令时是否可以进行任务切换?基于此 英特尔关于使用SSE3指令的文章 保存寄存器状态FSAVE和FRESTOR的命令已被FXSAVE和FXRESTOR取代,FXSAVE和FXRESTOR应该处理
累加器的全长。
在x64机器上,我认为“错误”指令可能包含在某些外部编译库中。
你当然会打 GCC Bug n°323,其他指出是由于FPU的精度过高。
方案:
- 使用SSE(或AVX,它是2016 ......)来执行计算
- 使用
-ffloat-store
编译开关。来自海湾合作委员会的文件。
不要将浮点变量存储在寄存器中,并禁止可能会改变从寄存器或内存中获取浮点值的其他选项。
这个选项可以防止诸如68000之类的机器上出现不希望的过度精度,其中浮动寄存器(68881)的精度比double应该具有的精度高。
同样适用于x86架构。
对于大多数程序来说,多余的精度只会很好,但有些程序依赖于IEEE浮点的精确定义。
在修改它们以将所有相关的中间计算存储到变量中之后,对这些程序使用-ffloat-store。