Taylor Series in C (problem with sin(240) and sin(300))

Question

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

const int TERMS = 7;
const float PI =  3.14159265358979;

int fact(int n) {
    return n<= 0 ? 1 : n * fact(n-1);
}


double sine(int x) {
    double rad = x * (PI / 180);
    double sin = 0;


    int n;
    for(n = 0; n < TERMS; n++) { // That's Taylor series!!
        sin += pow(-1, n) * pow(rad, (2 * n) + 1)/ fact((2 * n) + 1);
    }
    return sin;
}


double cosine(int x) {
    double rad = x * (PI / 180);
    double cos = 0;

    int n;

    for(n = 0; n < TERMS; n++) { // That's also Taylor series!
         cos += pow(-1, n) * pow(rad, 2 * n) / fact(2 * n);
    }
    return cos;
  }

int main(void){
   int y;
   scanf("%d",&y);
   printf("sine(%d)= %lf\n",y, sine(y));
   printf("cosine(%d)= %lf\n",y, cosine(y));

  return 0;
}

The code above was implemented to compute sine and cosine using Taylor series. I tried testing the code and it works fine for sine(120). I am getting wrong answers for sine(240) and sine(300).

Can anyone help me find out why those errors occur?

7 isn't probably enough for 240 or 300. But you can't go farther with integer factorial unless you use long long — Jean-François Fabre, Oct 23 '19 at 18:16
Note that [a POSIX-compliant system will define `M_PI` to the value of pi in ``](https://pubs.opengroup.org/onlinepubs/9699919799/basedefs/math.h.html#tag_13_26). — Andrew Henle, Oct 23 '19 at 18:23
You should forget about degree conversions when doing mathematical analysis, and work "natively" in radians. — Kaz, Oct 23 '19 at 20:10
`pow` can be very bad too... you want to calculate the powers by just multiplying in the loop. — Antti Haapala -- Слава Україні, Oct 23 '19 at 20:47

score 5 · Answer 1 · answered Oct 23 '19 at 20:21

You should calculate the functions in the first quadrant only [0, pi/2). Exploit the properties of the functions to get the values for other angles. For instance, for values of x between [pi/2, pi), sin(x) can be calculated by sin(pi - x).

The sine of 120 degrees, which is 40 past 90 degrees, is the same as 50 degrees: 40 degrees before 90. Sine starts at 0, then rises toward 1 at 90 degrees, and then falls again in a mirror image to zero at 180.

The negative sine values from pi to 2pi are just -sin(x - pi). I'd handle everything by this recursive definition:

sin(x):
  cases x of:
    [0, pi/2)   -> calculate (Taylor or whatever)
    [pi/2, pi)  -> sin(pi - x)
    [pi/2, 2pi) -> -sin(x - pi)
    < 0         -> sin(-x)
    >= 2pi      -> sin(fmod(x, 2pi))  // floating-point remainder

A similar approach for cos, using identity cases appropriate for it.

It is also a non-trivial issue that ther naive servies used in the OP's code is simply slow to converge. Chosing a faster converging series would help. — dmckee --- ex-moderator kitten, Oct 23 '19 at 22:04
@dmckee, that's the reason for the recommendation made in this answer. If you limit the absolute value of the angle below one radian, the convergence of the series is very quick, becoming slower as you depart from the origin value (in the case shown `0.0`) — Luis Colorado, Oct 26 '19 at 19:26

Jean-François Fabre · Answer 2 · 2019-10-23T19:28:58.037

The key point is:

TERMS is too small to have proper precision. And if you increase TERMS, you have to change fact implementation as it will likely overflow when working with int.

I would use a sign to toggle the -1 power instead of pow(-1,n) overkill.

Then use double for the value of PI to avoid losing too many decimals

Then for high values, you should increase the number of terms (this is the main issue). using long long for your factorial method or you get overflow. I set 10 and get proper results:

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

const int TERMS = 10;
const double PI =  3.14159265358979;

long long fact(int n) {
    return n<= 0 ? 1 : n * fact(n-1);
}
double powd(double x,int n) {
    return n<= 0 ? 1 : x * powd(x,n-1);
}


double sine(int x) {
    double rad = x * (PI / 180);
    double sin = 0;


    int n;
    int sign = 1;
    for(n = 0; n < TERMS; n++) { // That's Taylor series!!
        sin += sign  * powd(rad, (2 * n) + 1)/ fact((2 * n) + 1);
        sign = -sign;
    }
    return sin;
}


double cosine(int x) {
    double rad = x * (PI / 180);
    double cos = 0;

    int n;
    int sign = 1;
    for(n = 0; n < TERMS; n++) { // That's also Taylor series!
         cos += sign * powd(rad, 2 * n) / fact(2 * n);
         sign = -sign;
    }
    return cos;
  }

int main(void){
   int y;
   scanf("%d",&y);
   printf("sine(%d)= %lf\n",y, sine(y));
   printf("cosine(%d)= %lf\n",y, cosine(y));

  return 0;
}

result:

240
sine(240)= -0.866026
cosine(240)= -0.500001

Notes:

my recusive implementation of pow using successive multiplications is probably not needed, since we're dealing with floating point. It introduces accumulation error if n is big.
fact could be using floating point to allow bigger numbers and better precision. Actually I suggested long long but it would be better not to assume that the size will be enough. Better use standard type like int64_t for that.
fact and pow results could be pre-computed/hardcoded as well. This would save computation time.

Also, note some of the points brought up on this question and answers regarding greater floating-point precision: https://stackoverflow.com/questions/13516476/long-double-gcc-specific-and-float128 — Andrew Henle, Oct 23 '19 at 18:27
I think your changes are excessive (albeit likely improvements). The core problem with OP's code is the overflow in `fact`, and possibly also insufficient `TERMS`. — R.. GitHub STOP HELPING ICE, Oct 23 '19 at 19:24
Yes, you're right. My notes at the end of the answer explain that. The key point is : if TERMS is small then there's loss of precision, and if TERMS is >= 7 you get overflow unless you use int64_t or double. If you think you can edit the answer to make it more focused on the real issue, I see no issue with that. — Jean-François Fabre, Oct 23 '19 at 19:25
The coefficients could be precomputed in their entirety too, not just fact. — R.. GitHub STOP HELPING ICE, Oct 23 '19 at 19:27
true! an example where ideas of improvements came when writing the answer — Jean-François Fabre, Oct 23 '19 at 19:29
Thanks for the all the help. I have implemented my own power function and changed the the fact function to return a long long and these changes solved my issues. — FNF, Oct 24 '19 at 11:46

score 0 · Answer 3 · answered Oct 25 '19 at 12:31

const double TERMS = 14;
const double PI = 3.14159265358979;

double fact(double n) {return n <= 0.0 ? 1 : n * fact(n - 1);}

double sine(double x)
{
    double rad = x * (PI / 180);
    rad = fmod(rad, 2 * PI);
    double sin = 0;

    for (double n = 0; n < TERMS; n++) 
        sin += pow(-1, n) * pow(rad, (2 * n) + 1) / fact((2 * n) + 1);

    return sin;
}

double cosine(double x)
{
    double rad = x * (PI / 180);
    rad = fmod(rad,2*PI);
    double cos = 0;

    for (double n = 0; n < TERMS; n++)
        cos += pow(-1, n) * pow(rad, 2 * n) / fact(2 * n);

    return cos;
}

 int main()
{

    printf("sine(240)= %lf\n",  sine(240));
    printf("cosine(300)= %lf\n",cosine(300));
}

Taylor Series in C (problem with sin(240) and sin(300))

3 Answers3