C++实现Taylor Series计算cos(x) |x|>7时结果偏差问题

基于泰勒级数编写C++程序计算cos(x)值时，输入参数a在[-7,7]区间内运行结果正常，当a>7或a<-7时，计算结果和标准值偏差极大。
原问题代码如下：

#include <iostream>
#include <cmath>

using namespace std;


unsigned long long silnia(int n)
{
    unsigned long long a = 1;
    while(n)
        a *= n--;
    return a;
}
void taylor(int x)
{
    long double suma=0;
    for (int n=0; n<=10; n++)
    {
        suma+=pow(-1,n)*(pow(x,2*n)/silnia(2*n));
        cout << pow(-1,n) << " * " << pow(x,2*n) <<" / " << silnia(2*n) << " = " << pow(-1,n)*(pow(x,2*n)/silnia(2*n))<<endl;
    }
    cout << "taylor: "<<suma<< endl;
}

int main() {
    int a=5;
   taylor(a);
    cout << "cos: " << cos(a);
}

测试结果：

• a=5时，泰勒级数计算结果为0.283664，标准库cos函数计算结果为0.283662，误差符合预期：

taylor: 0.283664
cos: 0.283662

• a=9时，逐项计算输出如下，最终泰勒级数计算结果为-0.149111，标准库cos(9)正确结果为-0.91113，偏差明显：

1 * 1 / 1 = 1
-1 * 81 / 2 = -40.5
1 * 6561 / 24 = 273.375
-1 * 531441 / 720 = -738.112
1 * 4.30467e+007 / 40320 = 1067.63
-1 * 3.48678e+009 / 3628800 = -960.864
1 * 2.8243e+011 / 479001600 = 589.621
-1 * 2.28768e+013 / 87178291200 = -262.414
1 * 1.85302e+015 / 20922789888000 = 88.5647
-1 * 1.50095e+017 / 6402373705728000 = -23.4436
1 * 1.21577e+019 / 2432902008176640000 = 4.99719
taylor: -0.149111
cos: -0.91113

问题根因 #

代码存在两个核心问题，导致大输入下结果偏差：

1. 未对输入弧度做区间归约
   cos(x)是周期为2π的周期函数，且具备奇偶性、多段对称特性。当前使用的是在0点展开的麦克劳林级数（泰勒级数在0点的特殊形式），离展开点x=0越远，级数收敛越慢，达到相同精度需要计算的项数越多。代码固定只计算前11项（n取0到10），当|x|>7时已经离0点过远，11项远不足以收敛到正确值。另外9弧度已经超过2π（约6.28），完全可以利用周期性把值映射到离0点更近的区间，大幅降低需要计算的项数。
2. 单独计算高次幂与阶乘的方式存在数值溢出和精度损失
   每一项都单独计算pow(x,2n)和(2n)!：一方面unsigned long long的存储上限约为1.8e19，20!已经达到2.4e18，再往高阶计算阶乘会直接溢出，得到错误值；另一方面高次幂和大阶乘都是量级极大的数，两个大数做浮点除法会丢失大量有效精度，进一步放大误差。从a=9的逐项输出可以看到，算到n=10时项的绝对值还有4.99，远没有收敛到可以忽略的程度，此时截断累加自然误差极大。

修复方案 #

• 首先做输入值归约：利用cos函数的周期性、奇偶性，把任意输入的x先映射到[-π, π]甚至更小的[0, π/4]区间，让待计算值尽可能靠近0点，从根源上降低级数收敛需要的项数。比如x=9，减去2π后得到约2.717，cos(9)和cos(2.717)完全相等，后者离0点近很多，计算效率和精度都会大幅提升。
• 优化项的计算方式，放弃单独计算高次幂和阶乘：改用递推方式计算每一项，递推公式为第n项 = 第n-1项 * (-x*x) / ((2n-1)*2n)，n=0时初始项为1。这种方式全程不会产生极大的中间值，既避免了阶乘溢出问题，也大幅提升了浮点计算的数值稳定性，不需要额外写阶乘函数。
• 取消固定循环次数的写法：改为动态判断当前累加项的绝对值，当项的绝对值小于预设精度阈值（比如1e-10）时再停止累加，避免因为截断项数不足导致的误差。

修复后的参考代码：

#include <iostream>
#include <cmath>
using namespace std;

const long double PI = acosl(-1.0L);
const long double EPS = 1e-10L;

long double taylor_cos(long double x) {
    // 奇偶性归约：cos(-x)=cos(x)，统一转成非负值
    if (x < 0) x = -x;
    // 周期归约：cos(x+2kπ)=cos(x)，把x归约到[0, 2π)
    x = fmodl(x, 2*PI);
    // 进一步归约到[0, π]，利用cos(π+x)=-cos(x)
    int sign = 1;
    if (x > PI) {
        x -= PI;
        sign = -1;
    }
    long double suma = 0.0L;
    long double term = 1.0L; // n=0时的初始项
    suma += term;
    for (int n=1; ; n++) {
        // 递推计算当前项，无需单独计算pow和阶乘
        term = term * (-x*x) / ((2*n -1) * 2*n);
        suma += term;
        // 项足够小时退出循环，达到精度要求
        if (fabsl(term) < EPS) break;
    }
    return sign * suma;
}

int main() {
    int a = 9;
    long double res = taylor_cos(a);
    cout << "taylor: " << res << endl;
    cout << "cos: " << cos(a) << endl;
    return 0;
}

上述代码运行后，a=9时泰勒计算结果和标准库cos结果误差小于1e-10，即使输入更大的x值也能保持精度稳定。

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