C 中的泰勒级数 - math.h 与自己的实现

Posted 2023-03-28

【中文标题】C 中的泰勒级数 - math.h 与自己的实现

【英文标题】：Taylor series in C - math.h vs own implementation

【发布时间】：2021-11-13 16:41:36

【问题描述】：

虽然我在这里找到了一些有关泰勒系列的帖子，但我想寻求支持： 我根据泰勒方程为sin(x) 计算编写了C 代码。 作为参数，我的函数采用 rad 值和泰勒级数的预期长度。

我观察到的是，我的函数从math.h 返回与正弦相同的值，直到x <= 1.8 -> 所有高于此值的返回值都不同。代码如下：

这是在线调试器中的代码（这是我第一次粘贴，所以可能无法正常工作） https://onlinegdb.com/f4ymuMloW

完整代码：

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

int factorial(int x) 
    if (x == 0 || x == 1) 
        return 1;

    return x * factorial(x-1);


/*
    sin(x) = sum[k=0, n] ((-1^k) * x^(2k+1))/(2k+1)!        
 */

double taylorSine(double x, int k) 
    double ret = 0.0f;
    for (int i = 0 ; i < k; i++) 
        double l = pow(-1, i) * pow(x, 2 * i + 1);
        long long unsigned int m = factorial(2 * i + 1);
        ret += l/(double)m;
    
    return ret;


int main(int argc, char **argv) 
    float low,up;
    /* default */
    int sumsteps = 5;
    double res = 0.01f;

    if (argc == 1) 
        low = -3.14f;
        up = 3.14f;            
     
    else if (argc == 3) 
        low = atof(argv[1]);
        sumsteps = atoi(argv[2]);
     
    else 
        printf("wrong number of args\n");
        return 0;
    

    double r = taylorSine(low, sumsteps);
    printf("%f,%f\n", low, r);
    printf("sin(x)=%f\n", sin(low));        

    return 0;

和输出：

Starting program: /home/a.out 0 7
0.000000,0.000000
sin(x)=0.000000
[Inferior 1 (process 2146) exited normally]
(gdb) run 1.57 7
Starting program: /home/a.out 1.57 7
1.570000,1.000000
sin(x)=1.000000
[Inferior 1 (process 2147) exited normally]
(gdb) run 3.14 7
Starting program: /home/a.out 3.14 7
3.140000,0.002643
sin(x)=0.001593
[Inferior 1 (process 2148) exited normally]
(gdb) run 6.28 7
Starting program: /home/a.out 6.28 7
6.280000,9.053029
sin(x)=-0.003185
[Inferior 1 (process 2149) exited normally]
(gdb) 

【问题讨论】：

首先将factorial 改为使用unsigned long long 而不是int。

尝试更大的 sumsteps。泰勒级数的精度随着 x 的大小迅速恶化

sin 的良好实现不仅仅评估泰勒级数。他们将参数以 2π 为模（通常分为两部分，实际上是扇区数和以 π 的某个分数为模的残差），然后应用工程多项式，可能是极小极大多项式（选择以最小化相对误差而不是绝对误差）。在这些步骤中可能会使用一些扩展精度。

当您从上一次迭代中知道这些数字并且一些简单的数学（较少的计算机密集型）可以从前一个计算出下一个值时，为什么还要使用 pow 和 factorial！

@EdHeal 感谢您的建议，但是我从来不擅长这种优化，所以我坚持重复相同的操作，尽管我知道 x^7 = x^5 (prev) * x^2 和很快 ：） 。关于阶乘 - 我的想法只是创建查找表，但如果你有更快/更好/最佳的东西 - 你能给我一些提示吗？

【参考方案1】：

您的阶乘函数将溢出 (2*7+1)！ == 1,307,674,368,000，试试双版本吧。

double factorial(double x)

  double sum = 1.0;
  for (double y = 2.0; y < x; y += 1.0 )
  
    sum*=y;
  
  return sum;