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

Question

我试图在 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，如果它提供了任何有用的信息。

Answer 1

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

即使您使用了 long double为此，我已经在评论中指出 long double在 MSCRT 中映射到 double ( Reference ) 你仍然会在更大的角度看到一些异常，因为尽管 double可以存储与 1.8E+308 一样大的值但它在 2^53 = 9007199254740992 ~= 18! 处失去了粒度(即 2^53 + 1 存储为 double 等于 2^53 )。因此，一旦你的角度上升，这种行为的影响就会变得越来越大，以至于你在 printf() 中使用的小数点后 6 位精度很明显。 .

虽然您的方向是正确的，但您应该使用像 GMP 或 libcrypto 这样的 bignum 库。他们可以在不损失精度的情况下执行这些计算。

顺便说一句，由于您是在 Windows 7 上进行开发，这意味着您使用的是 x86 或 x86-64。在这些平台上，x87 能够使用 80 位执行扩展精度(根据 754 标准)操作，但我不知道编译器内部函数可以在不借助汇编代码的情况下为您提供该功能。

我还想提醒您注意范围缩小技术。虽然我仍然推荐使用 bignum 库，如果你在 0 之间很好的话和 90度(0 和 45 如果我要更严格的话)，你可以计算 sine()所有其他角度仅通过简单的三角恒等式。

更新:

实际上我要纠正自己关于使用 double 的问题s 在阶乘计算中。在编写了一个简单的程序后，我验证了当我使用 double 时存储阶乘，即使我高于 18，它们也是正确的.经过一番思考后，我意识到在阶乘的情况下，double 的情况。的粒度有点复杂。我举个例子来说明一下:

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

在这个号码18, 16, 14, ... , 2都是2的倍数自从乘以 2相当于在二进制表示中左移，19! 中的所有低位已经0因此当double对于大于 2^53 的 整数 开始舍入，这些阶乘不受影响。您可以计算 19! 的二进制表示中最低有效零的数量。通过计算 2 的数量这是16 . (对于 20! ，它是 18 )

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

更新 2:

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

对于 22! , 它是 51位(你可以自己计算)。