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);
}

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

Answer 1

每当for循环进行时，n增加2，因此对于dept=16，在接近循环结束时，必须计算大到30的数字的阶乘，并且使用的是unsigned int，它只能存储大到2^32=4294967296~=12！的值，这会导致阶乘函数溢出，从而给出错误的阶乘。

即使您对它使用了Long Double，并且我已经在注释中指出，MSCRT中的Long Double映射到Double(Reference)，您仍然会看到一些异常，可能是从更大的角度来看，因为尽管Double可以存储像1.8e+308，但它在2^53=9007199254740992~=18！中丢失了粒度（即2^53+1存储为因此，一旦增加角度，这种行为的影响就会变得越来越大，在printf()使用的6小数精度中可以明显看出。

尽管您已经走上了正确的轨道，但您应该使用像GMP或libcrypto这样的bignum库。他们可以在不损失精度的情况下执行这些计算。

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

如果我们使用double来保存阶乘，那么从23！开始，阶乘就会受到舍入的影响。它可以很容易地显示出来，因为2^74<23！<2^75，这意味着至少需要75位精度才能表示它，但是由于23！有19最低有效位，其值为0，因此它需要75-19=56，这比Double提供的53位要大。

对于22！，它是51位（您可以自己计算）。