equation solving - Why does my sin calculating code in C return the wrong value?

Question

I want to calculate the sinus of user inputs using my own functions based on this equation:

sin(x) = sum_(i=0)^n (-1)^i * (x^(2 i + 1)/((2 i + 1)!))

I have this code and to my understandings I do exactly the same as what's written in the equation:

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

int faculty(int factor)
{
    int result = 1;

    if (factor > 0)
    {
        for (int i = 1; i <= factor; i++)
        {
            result = result * i;
        }
    }
    else
    {
        result = 1;
    }

    return result;
}

double my_sin(double x, int n)
{
    double my_sin = 0;
    for (int j = 0; j < n ; j++)
    {
        double i = (double)j;
        double faculty_j = (double)faculty(2*j+1);
        my_sin = my_sin + (pow((-1.0), i) * (pow(x, (double)(2.0 * i + 1.0)) / faculty_j));
    }

    return my_sin;
}

int main()
{
    int n = 0;
    double x = 0;

    printf("x=");
    scanf("%lf", &x);

    printf("n=");
    scanf("%i", &n);

    printf("sin(%i)=", (int)x);
    printf("%lf
", my_sin(x, n));

    return 0;
}

However for example when I use x = 8 and n = 5 I get sin(8)=149 as a result. I tried debugging the code for some time now but I have no idea where the problem might be or how to find out what the problem is.

---

Updated code:

#include <stdio.h>

long int factorial(int factor)
{
    long int result = 1;

    if (factor > 0)
    {
        for (int i = 1; i <= factor; i++)
        {
            result = result * i;
        }
    }
    else
    {
        result = 1;
    }

    return result;
}

double my_pow(double a, double b)
{
    if (b == 0)
    {
        return 1;
    }

    double result = a;
    double increment = a;
    double i, j;

    for (i = 1; i < b; i++)
    {
        for (j = 1; j < a; j++)
        {
            result += increment;
        }
        increment = result;
    }
    return result;
}

double my_sin(double x, int n)
{
    double my_sin = 0;
    for (int j = 0; j < n ; j++)
    {
        double i = (double)j;
        double faculty_j = (double)factorial(2*i+1);
        my_sin = my_sin + (my_pow((-1.0), i) * (my_pow(x, 2.0 * i + 1.0) / faculty_j));
    }

    return my_sin;
}

int main()
{
    int n = 0;
    double x = 0;

    printf("x=");
    scanf_s("%lf", &x);

    printf("n=");
    scanf_s("%i", &n);

    printf("sin(%i)=", (int)x);
    printf("%lf
", my_sin(x, n));

    return 0;
}

question from:https://stackoverflow.com/questions/65870518/why-does-my-sin-calculating-code-in-c-return-the-wrong-value

Answer 1

The problem here is not necessarily your code. It's just your expectations on the equation. I tested the Taylor series for n=5 in Wolfram Alpha and subtracted sin(x) from the series, using the query

(sum (-1)^i * x^(2i+1)/(2i+1)!, i=0 to 5) - sin(x)  

to see the error.

https://www.wolframalpha.com/input/?i=%28sum+%28-1%29%5Ei+*+x%5E%282i%2B1%29%2F%282i%2B1%29%21%2C+i%3D0+to+5%29+-+sin%28x%29

As you can see, the series does a very good approximation when x is around 0, but when x goes above 3, the errors starts to be huge, and they are growing extremely fast. At x=8, the error has the size 67, which is obviously useless. To get a reasonable approximation, you need to use around n=15. Only five is way to little.

What you should do here is to take advantage of the fact that sin(x) = sin(x + k * 2 * PI) where k is an arbitrary integer. And if x is negative, you can just simple use the fact that sin(x) = -sin(-x). I'm showing a simple example of how it's done:

double my_sin(double x, int n)
{
    double my_sin = 0;
    double sign = 1.0;

    // sin(-x) = -sin(x)
    // This does not improve accuracy. It's just to make the rest simpler
    if(x<0) {
        x=-x;
        sign *= -1;
    }

    // sin(x) = sin(k*2*PI + x)
    x -= 2*PI*floor(x/(2*PI));

    // Continue as usual

    return sign * my_sin;
}

You want to be as close to x=0 as possible. In order to improve this even further, you can use the fact that sin(x) = -sin(x+PI), as shown below.

The above is enough to keep everything in the interval [0, PI], but we can actually do better than that by using some more clever math. We can use the fact that sin(x) is symmetric around x=PI/2 + k*PI. So sin(PI/2 -x) = sin(PI/2 + x). Utilizing that will keep you in the interval [0, PI/2] and when you're in that interval, n=2 is enough to get an error below 0.005. And for n=5, the error is below 10^-7. You can add this after the previous steps:

// sin(x) = -sin(x+PI)
if(x>PI) {
    x-=PI;
    sign *= -1;
}

// sin(PI/2 -x) = sin(PI/2 + x)
if(x>PI/2) 
    x = PI - x;

However, there is one bug in your code. Remember that the case j=n should be included, so change

for (int j = 0; j < n ; j++)

to

for (int j = 0; j <= n ; j++)

An alternative is to always call the function with 1 added to the argument.

And just for clarity. Remember that the Taylor series expects x to be in radians. If you're using degrees, you'll get the wrong result.

Just for completeness, here is the full function with a few other fixes that does not affect correctness but performance and readability:

double my_sin(double x, int n)
{
    double ret = 0;
    double sign = 1.0;

    if(x<0) {
        x=-x;
        sign *= -1;
    }

    x -= 2*PI*floor(x/(2*PI));

    if(x>PI) {
        x-=PI;
        sign *= -1;
    }

    if(x>PI/2) 
        x = PI - x;

    size_t denominator = 1;
    double numerator = x;
    int s = -1;

    for (int j = 0; j <= n ; j++) {
        denominator *= 2*j + 1;
        s *= -1;

        ret += s * numerator / denominator;

        numerator *= x*x;
    }

    return sign*ret;
}

Aiming for a certain precision

Say that you want 3 correct decimal digits. Then we can utilize the fact that the numerator grows slower than the denominator and do something like this:

    size_t denominator = 1;
    double numerator = x;
    double term;
    int s = -1;
    int j = 0;
    double tolerance = 0.001;

    do {
        denominator *= 2*j + 1;
        j++;
        s *= -1;

        term = numerator / denominator;

        ret += s * term;

        numerator *= x*x;
    } while(abs(term) > tolerance);