JosiP
JosiP

Reputation: 3079

Taylor series in C - math.h vs own implementation

Although I find some posts with Taylor series here, I'd like to ask for support: I wrote C code for sin(x) calculation, based on Taylor equation. As arguments my functions takes rad value and expected length of Taylor series.

What I'm observing is that my function returns same value as sine from math.h only till x <= 1.8 -> all returned values above this value differs. Here is the code:

Here is code in online debugger (it's first time I'm pasting, so it might be not working) https://onlinegdb.com/f4ymuMloW

whole code:

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

and outputs:

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) 

Upvotes: 0

Views: 766

Answers (1)

AndersK
AndersK

Reputation: 36102

Your factorial function will overflow (2*7+1)! == 1,307,674,368,000, try a double version instead.

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

Upvotes: 2

Related Questions