Visual C++ 2012 中的正弦计算不一致？

Posted 2023-02-17

【中文标题】Visual C++ 2012 中的正弦计算不一致？

【英文标题】：Sine computation inconsistency in Visual C++ 2012?

【发布时间】：2012-09-01 07:01:45

【问题描述】：

考虑以下代码：

// Filename fputest.cpp

#include <cmath>
#include <cstdio>

int main()

    double x;
    *(__int64 *) &x = 0xc01448ec3aaa278di64; // -5.0712136427263319
    double sine1 = sin(x);
    printf("%016llX\n", sine1);
    double sine2;
    __asm 
    fld x
    fsin
    fstp sine2
    
    printf("%016llX\n", sine2);
    return 0;

使用 Visual C++ 2012 (cl fputest.cpp) 编译并执行程序时，输出如下：

3FEDF640D8D36174
3FEDF640D8D36175

问题：

为什么这两个值不同？ 是否可以发布一些编译器选项以使计算的正弦值完全相同？

【问题讨论】：

一个操作可以完全在浮点寄存器中进行。当 80 位寄存器写入 64 位内存地址时，另一个会导致精度损失。 FSTP 文档说“将值存储在内存中时，会将值转换为单实数或双实数格式。”

fsin 方法使用具有 80 位精度的 x87 FPU，sin 在 MSVC（我使用的是 2010）中的实现似乎使用具有 128 位 xmm* 寄存器的 SSE . （另见this question。）

【参考方案1】：

这个问题不是由 long double 到 double 的转换引起的。这可能是由于数学库中的sin 例程不准确。

fsin 指令被指定为在其范围内的操作数（根据 Intel 64 和 IA-32 架构软件开发人员手册，2011 年 10 月，第 1 卷，第 8.3 卷）生成 1 ULP（长双精度格式）内的结果。 10），在四舍五入到最近的模式。在 Intel Core i7 上，fsin 的提问者的值，-5.07121364272633190495298549649305641651153564453125 或 -0x1.448ec3aaa278dp+2，产生 0xe.fb206c69b0ba402p-4。我们可以很容易地从这个十六进制看到最后 11 位是 100 0000 0010。这些是从 long double 转换时将被舍入的位。如果它们大于 100 0000 0000，则数字将向上取整。他们更大。因此，将这个 long double 值转换为 double 的结果是 0xe.fb206c69b0ba8p-4，等于 0x1.df640d8d36175p-1 和 0.93631021832247418590355891865328885614871978759765625。另请注意，即使结果低一个 ULP，最后 11 位仍将大于 100 0000 0000 并且仍会向上取整。因此，此结果在符合上述文档的 Intel CPU 上应该不会有所不同。

将此与直接计算双精度正弦进行比较，使用产生正确舍入结果的理想sin 例程。的值的正弦近似0.93631021832247413051857150785044253634581268961333520518023697738674775240815140702992025520721336793516756640679315765619707343171517531053811196321335899848286682535203710849065933755262347468763562（的Maple 10计算的）。最接近此的双精度为 0x1.df640d8d36175p-1。这与我们将 fsin 结果转换为 double 获得的值相同。

因此，差异不是由 long double 转换为 double 引起的；将 long double fsin 结果转换为 double 会产生与理想双精度 sin 例程完全相同的结果。

对于提问者的 Visual Studio 包使用的 sin 例程的准确性，我们没有规范。在商业图书馆中，允许 1 个 ULP 或多个 ULP 的错误是很常见的。观察正弦与双精度值舍入的点有多接近：距双精度值的距离为 0.498864 ULP（双精度 ULP），因此距舍入更改的点为 0.001136 ULP。因此，即使sin 例程中的一个非常轻微的错误也会导致它返回 0x1.df640d8d36174p-1 而不是更接近的 0x1.df640d8d36175p-1。

因此，我推测差异的根源是sin 例程中的一个非常小的错误。

【参考方案2】：

（注意：如 cmets 中所述，这在 VC2012 上不起作用。我把它留在这里作为一般信息。我不建议依赖任何依赖于优化级别的东西！）

我没有 VS2012，但是在 VS2010 编译器上，你可以在命令行上指定 /fp:fast，然后我得到相同的结果。这会导致编译器生成“快速”代码，这些代码不一定完全遵循 C++ 中所需的舍入和规则，但与您的汇编语言计算相匹配。

我不能在 VS2012 中尝试这个，但我想它有相同的选项。

这似乎也只能在优化的构建中使用/Ox 作为选项。

【参考方案3】：

见Why is cos(x) != cos(y) even though x == y?

作为David mentioned in the comment，差异来自将 FP 寄存器中的数据移动到不同大小的内存位置（寄存器/ram）。而且也不总是分配；即使是另一个附近的浮点操作也足以刷新 FP 寄存器，从而使任何试图保证特定值的尝试都徒劳无功。如果您需要进行比较，您可以通过将所有结果强制到内存位置来缓解其中的一些问题，如下所示：

float F1 = sin(a); float F2 = sin(b); if (F1 == F2)

但是，即使这样也可能行不通。最好的方法是接受任何浮点运算只会“大部分准确”，并且从程序员的角度来看，即使重复执行相同的运算，这个错误实际上也是不可预测的和随机的。而不是

if (F1 == F2)

你应该使用一些效果

if (isArbitrarilyClose(F1, F2))

或

if (absf(F1 - F2) <= n)

n 是一个很小的数字。

【讨论】：

My answer 证明这是不正确的：在这种情况下，将 long-double fsin 结果舍入为 double 不会产生不正确的值或与计算正确舍入的 @ 的结果不同的值987654330@直接。

浮点计算可能在某些甚至大多数情况下是一致的；但是，您必须非常小心，即使如此，编译器仍可能会阻碍您的努力。最好还是做好最坏的打算。

【参考方案4】：

VS2012 中的代码生成器进行了重大更改以支持automatic vectorization。部分更改是 x86 浮点数学现在在 SSE2 中完成，不再使用 FPU，这是必要的，因为 FPU 代码无法矢量化。 SSE2 使用 64 位精度而不是 80 位精度进行计算，因此很有可能由于四舍五入，结果相差一位。也是@J99在VS2010中可以和/fp:fast得到一致结果的原因，它的编译器仍然使用FPU，/fp:fast直接使用FSIN结果。

此功能非常有用，请查看链接网址中 Jim Hogg 的视频，了解如何利用它。

【讨论】：

64 位精度并不意味着 50% 的分布不正确。诸如此类的超越例程通常首先计算近似多项式的较低量值部分，最后添加最大分量，从而产生在很长一段时间内可以正确四舍五入的结果。