C语言手动实现sin泰勒级数与math.h库计算结果差异问题排查

• 截断误差问题：当前实现用的是正弦函数在0点的泰勒展开（麦克劳林展开），仅用有限项累加的话，拉格朗日余项会随输入x的绝对值增大而快速升高，项数固定时x越大误差越明显。测试用的7项展开仅在x绝对值小于2的区间能保证较高精度，x=6.28时余项已经远大于可接受阈值。
• 阶乘溢出问题：实现的factorial函数返回值为32位int类型，最大值仅为2^31-1≈2e9，而13!（对应7项展开的最高阶阶乘）已经达到6.2e9，超出int范围产生溢出，计算得到的阶乘值本身就是错误的，最终结果自然偏差极大。即使后续将阶乘转为unsigned long long存储，溢出已经在int计算阶段发生，无法修正。
• 计算精度损失问题：每一项单独调用pow计算x的高次幂和-1的幂次，通用幂函数本身存在额外精度损失，累加后误差进一步放大。

1. 输入参数范围规约 #

利用正弦函数的周期性、奇偶性和对称性，先将输入x缩放到[-π/2, π/2]区间再进行计算：

• 先用fmod(x, 2*M_PI)将x规约到[0, 2π)区间，消除周期带来的冗余
• 再通过对称性调整：若x>π则取x=x-2π，若x>π/2则取x=π-x并记录符号位，最终保证计算时x的绝对值不超过π/2≈1.57，刚好处于原有实现的高精度区间。

2. 改写项计算逻辑，避免阶乘溢出和pow调用 #

不要单独计算x的幂次和阶乘，通过递推的方式计算每一项的值：
第i项和第i+1项的关系为：term[i+1] = term[i] * (-x*x) / ((2*i+2)*(2*i+3))
初始项term[0] = x，累加时不需要调用任何幂函数、不需要单独计算阶乘，完全避免溢出问题，同时精度和计算效率都大幅提升。

3. 调整项数终止逻辑 #

可以放弃固定项数的输入，改为当当前累加项的绝对值小于1e-15（双精度精度阈值）时自动停止累加，既保证精度也不会做多余计算。

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

double taylorSine(double x, int max_steps) {
    // 第一步：参数规约到[-π/2, π/2]
    x = fmod(x, 2 * M_PI);
    int sign = 1;
    if (x > M_PI) {
        x -= 2 * M_PI;
    }
    if (x < -M_PI) {
        x += 2 * M_PI;
    }
    if (x > M_PI/2) {
        x = M_PI - x;
    }
    if (x < -M_PI/2) {
        x = -M_PI - x;
        sign = -1;
    }
    // 递推计算累加项
    double ret = x;
    double term = x;
    for (int i = 1; i < max_steps; i++) {
        term *= (-x * x) / ((2*i) * (2*i + 1));
        ret += term;
        // 可选：精度足够时提前退出
        if (fabs(term) < 1e-15) break;
    }
    return ret * sign;
}

int main(int argc, char **argv) {
    float low;
    int sumsteps = 7;

    if (argc == 3) {
        low = atof(argv[1]);
        sumsteps = atoi(argv[2]);
    } else if (argc != 1) {
        printf("wrong number of args\n");
        return 0;
    }

    double r = taylorSine(low, sumsteps);
    printf("自定义sin(%f) = %.10f\n", low, r);
    printf("标准库sin(%f) = %.10f\n", low, sin(low));        
    return 0;
}

修改后即使输入x=6.28，7项展开的结果也会和标准库几乎一致。

内容的提问来源于stack exchange，提问作者JosiP