sin，cos，tan不准确

为什么sinl在参数接近pi的非零倍数时会给出不正确的结果？为什么sinl在参数很大时会给出错误的结果？以下代码说明了这一点。sin，cos，tan不准确

请注意，用于初始化变量pi的数字并不完全匹配任何64位长的double值。编译器选择最接近的值，即3.14159265358979323851280895940618620443274267017841339111328125。预期的正弦值可以使用libquadmath，gnu MPFR lib或在线计算器（如http://www.ttmath.org/online_calculator）找到。

#include <stdio.h> 
#include <math.h> 

int main (int argc, char *argv []) 
    { 
    volatile long double pi = 3.14159265358979323846L; 
    volatile long double big = 9223372035086174241L; 
    volatile long double expected1 = -5.0165576126683320235E-20L; 
    volatile long double expected2 = -4.2053336735954077951E-10L; 
    double result; 
    double ex1 = expected1, ex2 = expected2; 

    result = sinl (pi); 
    printf("expected: %g, \nreturned: %g\n\n", ex1, result); 
    result = sinl (big); 
    printf("expected: %g, \nreturned: %g\n\n", ex2, result); 
    return 0; 
    }

我使用的是gcc 4.7.3。使用volatile会使编译器无法用硬编码结果替换sinl()调用。我的电脑有一个Intel Core i7处理器并运行Windows。我将结果打印为double而不是long double，因为我使用的gcc的mingw端口不支持打印long double。下面是程序输出：

expected: -5.01656e-020, 
returned: -5.42101e-020 

expected: -4.20533e-010, 
returned: -0.011874

来源

2013-06-02 ScottD

必需阅读：[每个计算机科学家应该了解的浮点运算_]（http://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html）（特别是针对您的第二个测试）。此外，它将有助于使用pi值精确到最后一位，而不是代码中的近似值。（不幸的是，'M_PI'是一个'double'，而不是'long double'，尝试使用'pi = acosl（-1.0L）;'）期望值来自哪里？ PI的正弦应该为0. –

@TedHopp：看在上帝的份上，如果你打算链接到那个页面并称之为“必需的阅读”，至少要阅读它。 – tmyklebu

@tmyklebu - 我读过它了。它出什么问题了？ –

这种不准确性可以追溯到由sinl库代码中使用的FSIN处理器指令。所述指令FSIN，FCOS，和FPTAN是不准确至1.0 ULP如Intel声称： http://notabs.org/fpuaccuracy/

来源

2013-06-02 06:29:04 ScottD

因此，如果我正确读取图形，至少部分的补救措施是使用模数来限制<0,2Pi）中的值，其中精度更高？ – Humungus

是的。但我认为以所需的精度进行模运算并不是直截了当的。一个挑战是一个长长的双倍不能很好地逼近pi。 – ScottD

GNU库的文档（可通过运行info libc math errors）列出了1个ULP “已知错误” 为cosl在x86和“x86_64的/ FPU “。它不记录sinl的任何内容。我可以在我的x86_64机器上以π/ 2为周围的cosl重现类似的巨大错误。

也许你应该将此作为一个文档错误报告给glibc和Linux手册页人员;我无法想象这是值得实施“正确的”修复。

如果你真的想要一个快速准确的sinl，我不太确定在哪里看。 CRlibm确实是sin（double的变体）。 MPFR会处理它，但它会比fsin慢很多倍。

来源

2013-06-02 08:07:48 tmyklebu

我给了libquadmath一个尝试。当然，在将其128位应答舍入到64位精度后，它很容易获得预期的结果。它比预期的处理器指令慢。我在pi附近的论点上测量的速度要慢15倍，而在大的论据上慢60倍。 – ScottD

@ScottD：很酷。感谢指针;我不知道libquadmath。 – tmyklebu

为了达到1个pi的倍数的ULP准确度，内部常数M_PI应该具有大约106位精度（或者对于长双倍为128位）。

在减少阶段，一个完美的实现将不得不以某种方式在减去(x - M_PI)后产生缺失的53或64位精度，因为一个幼稚的实现会将此中间值计算为零。当参数为非零的整数倍时，问题当然会变得越来越大。

对于1 ULP准确度，M_PI的内部精度为66位是不够的。然后再一次，我们可以重新阅读索赔，并检查1 ULP的准确性是否相对于结果或论证被索赔。

来源

2013-06-02 08:53:54

除非有一个聪明的方法来实现它，否则我认为你需要远远超过128位的pi来为trig函数进行参数缩减。计算具有64位精度的'0x1p16382' mod pi似乎需要pi到约16450位。也许再多一点，如果有一个整数倍的pi，特别好用'long double'来近似。 – tmyklebu

是的，一旦你超出了可以完美表现为双打的long long的范围，每个指数需要大约一个额外的比特。我只是从问题（x = 2^63-some）的背景中假设，人们并不真正期望得到准确的结果（10^6000）。 –

合理质量的三角函数的实现不使用'M_PI'来减少。使用Payne-Hanek算法。本质上，执行乘以1 /π，但它是一种专门的乘法，只返回小数位，丢弃整数值。由于整数值被丢弃，所以不需要1 /π的所有比特。所需的最高位由输入操作数的大小决定，最低位由一些数论性质决定（间隔中的某些数字与π/ 2的倍数有多接近）。 –

sin，cos，tan不准确

回答

相关问题