运行物理模拟大约20分钟后,错误陷阱将跳闸。意识到调试起来很麻烦,我在一个新项目中复制了相关子例程,并在发生错误时使用原始输入数据的硬编码副本对其进行了调用。但是错误陷阱没有跳闸!经过为期两天的繁琐工作,找出了子例程的两个实例的行为出现差异的确切点之后,我将问题追溯到了一个用于计算Cross_Product的非常简单的函数。
这两个程序的Cross_Product函数是相同的。我什至检查了反汇编,并确保编译器正在生成相同的代码。在两种情况下,该功能还接收相同的输入数据。我什至明确地检查了函数内部的舍入模式,它们是相同的。但是,他们返回的结果略有不同。具体来说,对于返回的三个 vector 分量中的两个,LSB是不同的。甚至调试器本身也确认三个变量中的这两个不等于它们已明确分配的表达式。 (请参见下面的屏幕截图。)
在原始程序中,调试器在监视列表的最后三行中均显示“true”,而不仅仅是最后一行。
我在XP和GCC编译器以及AMD Athlon 64 CPU上使用Code::Blocks 13.12。但是,我在具有Intel Core i5 CPU的更现代的Windows 10计算机上重新编译并在Code::Blocks 16.01中运行了测试程序,结果是相同的。
这是我最少,完整且可验证的代码,可再现与我的原始程序以及调试器本身不符的奇怪结果(不幸的是,由于它的庞大,我无法包含原始物理程序):
extern "C" {
__declspec(dllimport) int __stdcall IsDebuggerPresent(void);
__declspec(dllimport) void __stdcall DebugBreak(void);
}
struct POLY_Triplet {
double XYZ[3];
};
POLY_Triplet Cross_Product(POLY_Triplet Vector1, POLY_Triplet Vector2) {
POLY_Triplet Result;
Result.XYZ[0] = Vector1.XYZ[1] * Vector2.XYZ[2] - Vector1.XYZ[2] * Vector2.XYZ[1];
Result.XYZ[1] = Vector1.XYZ[2] * Vector2.XYZ[0] - Vector1.XYZ[0] * Vector2.XYZ[2];
Result.XYZ[2] = Vector1.XYZ[0] * Vector2.XYZ[1] - Vector1.XYZ[1] * Vector2.XYZ[0];
return Result;
}
int main() {
POLY_Triplet Triplet1;
POLY_Triplet Collision_Axis_Vector;
POLY_Triplet Boundary_normal;
*(long long int *)(&Collision_Axis_Vector.XYZ[0]) = 4594681439063077250;
*(long long int *)(&Collision_Axis_Vector.XYZ[1]) = 4603161398996347097;
*(long long int *)(&Collision_Axis_Vector.XYZ[2]) = 4605548671330989714;
*(long long int *)(&Triplet1.XYZ[0]) = -4626277815076045984;
*(long long int *)(&Triplet1.XYZ[1]) = -4637257536736295424;
*(long long int *)(&Triplet1.XYZ[2]) = 4589609575355367200;
if (IsDebuggerPresent()) {
DebugBreak();
}
Boundary_normal = Cross_Product(Collision_Axis_Vector, Triplet1);
return 0;
}
为了方便起见,以下是监视列表的相关行,如屏幕截图所示:
(Result.XYZ[0] == Vector1.XYZ[1] * Vector2.XYZ[2] - Vector1.XYZ[2] * Vector2.XYZ[1])
(Result.XYZ[1] == Vector1.XYZ[2] * Vector2.XYZ[0] - Vector1.XYZ[0] * Vector2.XYZ[2])
(Result.XYZ[2] == Vector1.XYZ[0] * Vector2.XYZ[1] - Vector1.XYZ[1] * Vector2.XYZ[0])
谁能解释这个行为?
最佳答案
我可以确认您得到的输出有问题是由x87精度的变化引起的。精度值存储在x87 FPU控制寄存器中,更改后,该值将在线程的整个生命周期中持续存在,从而影响线程上运行的所有x87代码。
显然,您的大型程序(或所使用的外部库)的其他某些组件有时会将尾数长度从53位(默认值)更改为64位(这意味着使用这80位x87寄存器的全精度)。
修复的最佳方法是将编译器从x87切换到SSE2目标。 SSE始终使用32位或64位浮点数(取决于所使用的指令),它根本没有80位寄存器。甚至您的2003 Athlon 64也已经支持该指令集。副作用是,您的代码将变得更快。
更新:如果您不想切换到SSE2,则可以将精度重置为任意值。在Visual C++中执行此操作的方法如下:
#include <float.h>
uint32_t prev;
_controlfp_s( &prev, _PC_53, _MCW_PC ); // or _PC_64 for 80-bit
对于海湾合作委员会来说,就是这样(未经测试)
#include <fpu_control.h>
#define _FPU_PRECISION ( _FPU_SINGLE | _FPU_DOUBLE | _FPU_EXTENDED )
fpu_control_t prev, curr;
_FPU_GETCW( prev );
curr = ( prev & ~_FPU_PRECISION ) | _FPU_DOUBLE; // or _FPU_EXTENDED for 80 bit
_FPU_SETCW( curr );