2012-09-01 18 views
8

考虑下面的代码:Visual C++ 2012中的正弦计算不一致性?

// 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 

问题:

  • 为什么这两个值不同?
  • 是否可以发出一些编译器选项,以便计算出的正弦值完全相同?
+1

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

+1

'fsin'方法使用80位精度的x87 FPU,在MSVC(我使用2010)中执行'sin'似乎使用SSE及其128位xmm *寄存器。 (另请参阅[此问题](http://stackoverflow.com/questions/2284860/how-does-c-compute-sin-and-other-math-functions)。) – DCoder

回答

9

此问题不是由long double转换为double导致的。这可能是由于数学库中的sin例程不准确。

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

比较此直接计算双精度正弦,使用理想的sin例程,产生正确的舍入结果。值的正弦值大约是0.93631021832247413051857150785044253634581268961333520518023697738674775240815140702992025520721336793516756640679315765619707343171517531053811196321335899848286682535203710849065933755262347468763562(用枫10计算)。与此最接近的是0x1.df640d8d36175p-1。这与我们通过将fsin结果转换为两倍而获得的值相同。

因此,该差异不是由长双倍转换为双倍所引起的;将长双倍结果转换为双倍结果产生与理想双精度例程完全相同的结果。

我们没有关于提问者的Visual Studio包使用的sin例程的准确性的规范。在商业图书馆中,允许1 ULP或几个ULP的错误是常见的。观察正弦到双精度值四舍五入的点的距离:从双精度值开始,它是.498864 ULP(双精度ULP),因此它距离舍入更改的点为0.001136 ULP。因此,即使sin例程中的一个非常微小的不准确性,也会导致它返回0x1.df640d8d36174p-1而不是更接近0x1.df640d8d36175p-1。

因此,我猜想这种差异的根源在sin例程中是一个非常小的不准确性。

1

(注意:!正如评论所说,这并不对VC2012工作,我在这里一般信息离开它我会建议不要依赖任何依赖于优化级别反正)

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

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

这似乎只在/Ox作为一个选项优化构建工作。

1

参见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是一个很小的数量的影响。

+2

[我的回答](http:/ /stackoverflow.com/a/12227672/298225)证明这是不正确的:在这种情况下,将long-double'fsin'结果舍入为double不会产生不正确的值或不同于计算正确舍入的'罪'直接。 –

+0

浮点计算*可能在某些甚至大多数情况下是一致的;然而,你必须非常小心,即使如此,编译器仍然可能会阻止你的努力。为最坏的情况做好准备仍然是最好的。 – Ghost2

0

VS2012中的代码生成器已大幅更改为支持automatic vectorization。这一变化的部分原因是x86浮点数学现在已经在SSE2中完成,不再使用FPU,因为FPU代码不能被矢量化。 SSE2以64位精度而不是80位精度计算,结果由于四舍五入导致结果偏离一位。也就是@ J99可以在VS2010中以/ fp:fast得到一致的结果,其编译器仍然使用FPU和/ fp:fast直接使用FSIN结果。

这个功能有很多功能,请查看Jim Hogg在链接网址上的视频,了解如何利用它。

+2

64位精度并不意味着错误结果的50%分布。先验例程,如这通常计算的近似多项式的第一低幅度部分和添加的最大组成部分最后,产生可能是正确的舍入的大量时间的结果。 –