C 中正弦函数的实现不起作用

我尝试用 C 语言实现正弦函数，但得到了奇怪的结果。以下是我用来计算正弦的三个函数：

#define PI 3.14159265358979323846
#define DEPTH 16

double sine(long double);
long double pow(long double, unsigned int);
unsigned int fact(unsigned int);

double sine(long double x) {
    long double i_x = x *= PI/180;
    int n = 3, d = 0, sign = -1;

    // fails past 67 degrees
    for (; d < DEPTH; n += 2, d++, sign *= -1) {
        x += pow(i_x, n) / fact(n) * sign;
    }

    return x;
}

long double pow(long double base, unsigned int exp) {
    double answer = 1;
    while (exp) {
        answer *= base;
        exp--;
    }
    return answer;
}

unsigned int fact(unsigned int n) {
    unsigned int answer = 1;
    while (n > 1) {
        answer *= n--;
    }
    return answer;
}

为了测试它，我已经针对内置正弦函数对其进行了测试，如下所示：

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

main() {
    for (int i = 0; i <= 180; i++) {
        printf("sin(%i) = %lf, %lf\n", i, sine(i), sin(i*3.14159265358979323846/180));
    }

    exit(EXIT_SUCCESS);
}

直至 67 度，其计算结果与内置函数相同。不过，当它超过 67 时，它通常会偏离实际值越来越远。

这是一个示例输出：

>> sin(100) = 0.987711, 0.984808
>> sin(101) = 0.986885, 0.981627
>> sin(102) = 0.987056, 0.978148
>> sin(103) = 0.988830, 0.974370
>> sin(104) = 0.993060, 0.970296
>> sin(105) = 1.000948, 0.965926
>> sin(106) = 1.014169, 0.961262
>> sin(107) = 1.035052, 0.956305
>> sin(108) = 1.066807, 0.951057
>> sin(109) = 1.113846, 0.945519
>> sin(110) = 1.182194, 0.939693
>> sin(111) = 1.280047, 0.933580
>> sin(112) = 1.418502, 0.927184
>> sin(113) = 1.612527, 0.920505
>> sin(114) = 1.882224, 0.913545
>> sin(115) = 2.254492, 0.906308
>> sin(116) = 2.765192, 0.898794
>> sin(117) = 3.461969, 0.891007
              ...
>> sin(180) = 8431648.192239, 0.000000

有谁知道为什么会发生这种情况？我正在 Windows 7 上使用 Visual Studio 2017，如果这提供了任何有用的信息。

每次你的for loop进步，n增加了2因此对于DEPTH = 16，在循环结束时，您必须计算大到的数字的阶乘30你正在使用unsigned int只能存储大小为2^32 = 4294967296 ~= 12!这会导致阶乘函数溢出，从而给出错误的阶乘。

即使你用过long double为此，我已经在评论中指出long double在 MSCRT 中映射到double (参考 https://msdn.microsoft.com/en-us/library/9cx8xs15.aspx）你仍然会在更大的角度看到一些异常，因为尽管double可以存储大到的值1.8E+308但它失去了它的粒度2^53 = 9007199254740992 ~= 18! (i.e. 2^53 + 1存储为double等于2^53）。因此，一旦角度上升，这种行为的影响就会变得越来越大，以至于在您使用的 6 位小数精度中很明显printf().

虽然你走在正确的轨道上，但你应该使用像这样的bignum库GMP or 加密库。他们可以在不损失精度的情况下执行这些计算。

顺便说一句，由于您正在 Windows 7 上进行开发，这意味着您要么使用 x86，要么使用 x86-64。在这些平台上，x87 能够执行 80 位的扩展精度（根据 754 标准）操作，但我不知道编译器内在函数可以为您提供这种功能，而无需求助于汇编代码。

我还想请您注意范围缩小技术。虽然我仍然建议使用 bignum libs，如果你擅长两者之间0 and 90度（0 and 45如果我要更严格的话），你可以计算sine()所有其他角度都可以通过简单的三角恒等式.

UPDATE:

实际上我会纠正自己的使用double阶乘计算中的 s。编写了一个简单的程序后，我验证了当我使用double存储阶乘，即使我高于18。经过一番思考后，我意识到在这种情况下阶乘，情况与double的粒度稍微复杂一些。我举个例子让你明白：

19! = 19 * 18 * ... * 2 * 1

在这个数字中18, 16, 14, ... , 2都是的倍数2并且由于乘以2相当于二进制表示中的左移，所有低位都在19!已经0因此当double的舍入开始为integers比...更棒2^53，这些阶乘不受影响。您可以计算二进制表示中最低有效零的数量19!通过计算2是16. (for 20!, it is 18)

我要上去1.8e+308并检查所有阶乘是否不受影响。我会向您通报最新结果。

更新2：

如果我们使用doubles 保存阶乘，它们受到四舍五入的影响23!向前。它可以很容易地显示出来，因为2^74 < 23! < 2^75这意味着至少需要 75 位精度来表示它，但由于23! has 19最低有效位的值为0，所以需要75 - 19 = 56它大于53位提供者double.

For 22!, it is 51位（您可以自己计算）。