What is a more accurate algorithm I can use to calculate the sine of a number?

Question

I have this code that calculates a guess for sine and compares it to the standard C library's (glibc's in my case) result:

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

double double_sin(double a)
{
    a -= (a*a*a)/6;

    return a;
}

int main(void)
{
    double clib_sin = sin(.13),
             my_sin = double_sin(.13);
    printf("%.16f\n%.16f\n%.16f\n", clib_sin, my_sin, clib_sin-my_sin);
    return 0;
}

The accuracy for double_sin is poor (about 5-6 digits). Here's my output:

0.1296341426196949
0.1296338333333333
0.0000003092863615

As you can see, after .12963, the results differ.

Some notes:

• I don't think the Taylor series will work for this specific situation, the factorials required for greater accuracy aren't able to be stored inside an unsigned long long.

• Lookup tables are not an option, they take up too much space and generally don't provide any information on how to calculate the result.

• If you use magic numbers, please explain them (although I would prefer if they were not used).

• I would greatly prefer an algorithm is easily understandable and able to be used as a reference over one that is not.

• The result does not have to be perfectly accurate. A minimum would be the requirements of IEEE 754, C, and/or POSIX.

• I'm using the IEEE-754 double format, which can be relied on.

• The range supported needs to be at least from -2*M_PI to 2*M_PI. It would be nice if range reduction were included.

What is a more accurate algorithm I can use to calculate the sine of a number?

I had an idea about something similar to Newton-Raphson, but for calculating sine instead. However, I couldn't find anything on it and am ruling this possibility out.

Answer 1

You can actually get pretty close with the Taylor series. The trick is not to calculate the full factorial on each iteration.

The Taylor series looks like this:

sin(x) = x^1/1! - x^3/3! + x^5/5! - x^7/7!

Looking at the terms, you calculate the next term by multiplying the numerator by x^2, multiplying the denominator by the next two numbers in the factorial, and switching the sign. Then you stop when adding the next term doesn't change the result.

So you could code it like this:

double double_sin(double x)
{
    double result = 0;
    double factor = x;
    int i;

    for (i=2; result+factor!=result; i+=2) {
        result += factor;
        factor *= -(x*x)/(i*(i+1));
    }
    return result;
}

My output:

0.1296341426196949
0.1296341426196949
-0.0000000000000000

EDIT:

The accuracy can be increased further if the terms are added in the reverse direction, however this means computing a fixed number of terms:

#define FACTORS 30

double double_sin(double x)
{
    double result = 0;
    double factor = x;
    int i, j;
    double factors[FACTORS];

    for (i=2, j=0; j<FACTORS; i+=2, j++) {
        factors[j] = factor;
        factor *= -(x*x)/(i*(i+1));
    }
    for (j=FACTORS-1;j>=0;j--) {
        result += factors[j];
    }
    return result;
}

This implementation loses accuracy if x falls outside the range of 0 to 2*PI. This can be fixed by calling x = fmod(x, 2*M_PI); at the start of the function to normalize the value.