6

我的基本问题是如何使 x86 上的浮点运算表现得像 PowerPC，从 Classic MacOS (CodeWarrior) 到 Windows (VS 2008)。

有问题的代码有很多，有一堆算法，这些算法高度迭代并且对数值非常敏感。

典型的复杂线是：

Ims_sd = sqrt((4.0*Ams*sqr(nz)-8.0*(Ams+Dms)*nz+12.0*sqr(Ams)) /
         (4.0*sqr(Ams)*(sqr(nz)-1)) - 
         sqr(Ims_av))*sqrt(nz-1);

它是使用 typedef'dfloat作为基本类型编写的。

更改为double在两个平台上得到非常相似的结果，但不幸的是这些数字是不可接受的，所以我们不能采取那么简单的方法。

Mac 代码是使用 CodeWarrior 编译的，只是关闭 FMADD 和 FMSUB 指令的生成对创建的数字产生了巨大影响。因此，我的出发点是搜索看起来最相似的 Visual Studio (2008) 选项 - 确保使用了 fused add。我们怀疑关键在于编译器在计算中分配中间存储的行为

目前，通过启用 SSE2 和/fp:fast. 启用内在函数会导致值偏离 Mac 值。

/ fp开关文档说只/fp:strict关闭融合添加行为。

MSDN谈到“在 LIBC.LIB、LIBCMT.LIB 或 MSVCRT.LIB 之前”链接 FP10.OBJ。保证64位精度。我显然已经通过在链接器输入字段上指定 FP10.OBJ 来实现这一点（详细的链接器输出在 MSVCRTD.lib 之前显示它）。

我还通过调用设置了 64 位精度

_controlfp_s(&control_word, _PC_64, MCW_PC);

在 DllMain 中。

请注意，问题不是由于平台之间浮点异常处理的差异，也不是由于 PowerPC 允许除以零整数（仅返回零）的（令人愉快的）方式，因为这些区域已经过审计和解决，非常感谢PC-皮棉。该程序运行并产生了一些看似合理的输出，但还不够好。

更新：

一位朋友的有趣评论： 一种可能是 PPC 有大量临时寄存器，可以存储 64 位中间值，而 x86 代码可能必须卸载和重新加载 FPU（截断到 4 个字节并丢失精度）。

这可能是 SSE2 工作得更好的原因，因为 (IIRC) 它有更多的寄存器和更多的保留中间值的空间。

一种可能性 - 您的代码可以编译为 64 位吗？x64 模式还有更多的中间寄存器，以及更好的 FP 指令，因此在设计和执行上可能更接近 PPC。

正如他所建议的那样，使用 64 位构建的初始测试实际上更接近了（我最初认为它过头了，但这是由于建模设置不正确造成的）。

最终决议

我敢肯定，任何对这个话题感兴趣的人都足够痴迷，他们想知道这一切最终是如何解决的。该软件已完成并提供一致的数字结果。我们永远无法获得所有算法来为 Mac 提供相同的结果，但它们足够接近，可以在统计上接受。鉴于处理是由专家用户选择感兴趣的区域指导的，并且用户输入对模型的进展有部分反应，首席科学家认为这是可以接受的（这不是一夜之间的决定！）。剩余的数字差异完全在决定不同临床结果的范围内，因此在测试中没有看到不同的诊断。

4

3 回答 3

3

跨多个平台的浮点确定性的整个问题似乎是一个非常棘手的问题，你越深入它，它似乎变得越糟。

我确实发现了这篇有趣的文章，它深入讨论了这个问题——它可能会提出一些想法。

于 2010-02-28T14:31:14.190 回答
1

不是这样的答案，而是比我可以在评论中容纳的更多文本（和格式）。读到你的问题，我觉得你可能已经考虑了所有这些，但没有告诉我们，所以这可能都是无关紧要的闲聊。如果是，我道歉。

你能（你是吗？）在程序的原始版本或移植版本上强制遵守 IEEE754 浮点运算规则？我的第一个猜测是这两个平台（硬件、o/s、库的组合）实现了 fp 算法的不同方法。

您对两个平台上某些基本类型（如整数和浮点数）的默认大小做出了哪些假设（如果有的话）。C 标准（我相信 C++ 标准）允许平台依赖于某些此类（我无法记住哪个，我真的是 Fortran 程序员）。

最后的猜测——我已经习惯（在我的 Fortranny 世界中）指定浮点常量，例如你的 4.0，它有足够的数字来指定首选表示中的所有（十进制）数字，例如 4.000000000000000000000000。我知道，在 Fortran 中，像 3.14159625 这样的 4 字节浮点常量在自动转换为 8 字节时，不会用 pi 的十进制表达式中的其他数字填充额外的字节。这可能会影响您。

这些都不能真正帮助您确保代码的移植版本产生与原始版本相同的结果，只是识别差异的来源。

最后，您是否要求新版本产生与旧版本相同的结果，或者您向您的客户保证新版本产生准确的答案？考虑到数值计算中的所有错误来源，您的问题留下了旧版本程序比新版本“错误”的可能性。

于 2010-01-28T09:31:34.360 回答
1

我向您推荐 GCC 错误 323：

我想欢迎 bug 323 社区的最新成员，在这里 gcc 中的所有 x87 浮点错误都会消失！所有使用 x87 的浮点错误都是受欢迎的，尽管它们中的许多很容易修复，而且很多都不是！我们都是一个幸福的家庭，犯了一个严重的错误，即从市场上最准确的通用 FPU 中获得准确性！

简短的总结是，在 x87 上获得“真正的”IEEE 浮点单数/双精度而没有显着的性能损失是非常乏味的。即使fldcw由于减小的指数范围（IIRC，IEEE FP 特别允许实现自己做 WRT denorms），您也会遭受 denorms 的双舍入。大概你可以做这样的事情：

  1. 四舍五入到正无穷，执行操作（得到 ldresult1），四舍五入到最接近的偶数，转换为浮点数（得到 fresult1）。
  2. RTNI，执行op，RTNE，转换为float。
  3. 如果它们相同，那就太好了：您有正确的 RTNE 浮点结果。如果不是，那么（我认为）fresult2 < fresult1，进而，fresult1=nextafterf(fresult2,+inf)，有两种可能：
    • ldresult1 == ((long double)fresult1+fresult2)/2。“正确”的答案是fresult2。
    • ldresult2 == ((long double)fresult1+fresult2)/2。“正确”的答案是fresult1。

我可能在某处的细节上错了，但这大概是当你得到一个 denorm 时你必须经历的痛苦。

然后你遇到了另一个问题：我很确定 sqrt() 不能保证在不同的实现中返回相同的结果（并且非常确定对于 trig 函数）；我见过的唯一保证是结果“在 1 ulp 以内”（大概是正确舍入的结果）。它高度依赖于所使用的算法，并且现代 CPU 有这些指令，因此如果您尝试在软件中实现它，您将遭受显着的性能损失。尽管如此，ISTR 是一个“可移植的”浮点库，它应该可以实现一致性，但我不记得 OTTOMH 的名字了。

于 2010-10-03T03:15:57.207 回答