Calculate the function sin()

Question

For my studies, I have to code an algorithm to calculate sin() with this function:

However, in my algorithm, I have to keep the value of X between 0 and Pi/2. So, I wrote my algorithm but all the results are wrong.

Here is my code:

double sinX(double x){
    double resultat = 0;
    int i;
    if(x < 0 || x > M_PI_2)
        x = fmod(x,M_PI_2);

    for(i = 1;i<=30;i++){
        resultat += -1 * ((x*x)/(2*i*(2*i+1)))*(pow(-1,i-1))*((pow(x,2*i-1))/(factorielle(2*i-1)));
    }

    return resultat;
}

I didn't find the reason. Can you help me?

Here the are few values of X and the result with fmod

1 / 1
2 / 0.429204
3 / 1.4292
4 / 0.858407
5 / 0.287611
6 / 1.28761
7 / 0.716815
8 / 0.146018
9 / 1.14602
10 / 0.575222
11 / 0.00442571
12 / 1.00443
13 / 0.433629
14 / 1.43363
15 / 0.862833
16 / 0.292037
17 / 1.29204
18 / 0.72124
19 / 0.150444
20 / 1.15044

and the result with the algorithm

1 / -0.158529
2 / -0.0130568
3 / -0.439211
4 / -0.101605
5 / -0.00394883
6 / -0.327441
7 / -0.0598281
8 / -0.000518332
9 / -0.234888
10 / -0.0312009
11 / -1.44477e-008
12 / -0.160572
13 / -0.0134623
14 / -0.443022
15 / -0.103145
16 / -0.00413342
17 / -0.330639
18 / -0.0609237
19 / -0.000566869
20 / -0.237499

Here is my "factorielle" definition

double factorielle(double x){
    double resultat = 1;
    int i;

    if(x != 0){
        for (i=2;i<=x;i++)
        {
            resultat *= i;
        }
    }
    else{
        resultat = 1;
    }

    return resultat;
}

And values :

1 / 1
2 / 2
3 / 6
4 / 24
5 / 120
6 / 720
7 / 5040
8 / 40320
9 / 362880
10 / 3.6288e+006
11 / 3.99168e+007
12 / 4.79002e+008
13 / 6.22702e+009
14 / 8.71783e+010
15 / 1.30767e+012
16 / 2.09228e+013
17 / 3.55687e+014
18 / 6.40237e+015
19 / 1.21645e+017
20 / 2.4329e+018

Answer 1

You're misunderstanding the purpose of the second formula you show. The idea is that you use that formula to compute each term in the sum from the preceding term, saving you from the need to use any pow or factorial calls.

#include <stdio.h>

double sinX(double x) {
  double term, total_so_far;
  int i;

  term = x;  /* First term in the expansion. */
  total_so_far = 0.0;
  for (i = 1; i <= 30; i++) {
    /* Add current term to sum. */
    total_so_far += term;
    /* Compute next term from the current one. */
    term *= -(x * x) / (2*i) / (2*i + 1);
  }
  return total_so_far;
}


int main(void) {
  /* testing */
  double x;
  int i;

  for (i = 0; i <= 10; i++) {
    x = i / 10.0;
    printf("sin(%f) is %f\n", x, sinX(x));
  }
  return 0;
}

And the results of running this code, on my machine:

sin(0.000000) is 0.000000
sin(0.100000) is 0.099833
sin(0.200000) is 0.198669
sin(0.300000) is 0.295520
sin(0.400000) is 0.389418
sin(0.500000) is 0.479426
sin(0.600000) is 0.564642
sin(0.700000) is 0.644218
sin(0.800000) is 0.717356
sin(0.900000) is 0.783327
sin(1.000000) is 0.841471

That should give you reasonable results for the range 0 to pi / 2. Outside that range you'll need to be a bit cleverer about the reduction you're using: simply reducing modulo pi / 2 won't give correct results. (Hint: it's safe to reduce modulo 2 * pi, since the sin function is periodic with period 2 * pi. Now use symmetries of the sin function to reduce to the range 0 to pi / 2.)

---

EDIT An explanation of why the current code is giving incorrect results: apart from the flawed reduction step, in your sum you start with the term i = 1. But the first term should be for i = 0 (that's the x term, while the i=1 term is the -x^3 / 3! term). A quick and dirty fix is to remove the reduction step, and to initialise your resultat variable to x rather than 0. That should give you good results for small x, and then you can figure out how to replace the reduction step. I'd be surprised if you were really intended to compute the answer using explicit factorial and power calls, though - I'm almost sure that you're expected to compute each term from the preceding one as described above.

Answer 2

There are two problems with your code:

1. sin(x+k*π/2) does not necessarily equal sin(x)
2. Your expression for the term is a little bit messed up. The instructions seem to suggest that you calculate the next term in the series from the previous term. Start with the value for i=0 and then use the equation in your question to compute the next term in each iteration.

Answer 3

Finally, I followed yours indications. Here is my final code :

double sinX(double x)
{
    double result = 1.0;
    double term_i = 1.0;
    int i = 2;

    x = fmod(x, 2*M_PI);

    for(i = 2; i<= 30; i+=2)
    {
        term_i = (-term_i * (x*x)) / (i*(i+1));
        result += term_i;
    }

    return x * result;
}

Answer 4

Idea about the number of terms with OP's posted answer.

As long as one performs some range limitation first, like fmod(), the number of terms needed can be reasonably determined dynamically. (Uses 1 to 23 iterations for x: 0 to 2*pi.)

double sinX1(double x)
{
    double result = 1.0;
    double term_i = 1.0;
    int i = 2;

    x = fmod(x, 2*M_PI);

    // for(i = 2; i<= 30; i+=2)
    for(i = 2; ((1.0 + term_i) != 1.0); i+=2)
    {
        term_i = (-term_i * (x*x)) / (i*(i+1));
        result += term_i;
    }
    return x * result;
}