Incremental Sine/Cosine Generation

Steve Hollasch

I've coded a great many loops that sweep from 0 to 360 degrees. They all follow this basic pattern:

double theta_increment = (2*PI) / nslices;
int    i;
double Tcos, Tsin;

for (i=0;  i < nslices;  ++i)
{
    double angle = i * theta_increment;

    Tcos = cos (angle);
    Tsin = sin (angle);

    ...     // Use Tcos and Tsin
}

This loop generates the sine and cosine of an incremental angle. I might, for example, be spinning out vertices for a surface of revolution. Or I might be incrementally evaluating some function of sine or cosine. The question I had was whether I could improve such a loop to take advantage of the fact that I am using a constant increment, and possibly even eliminate the calls to sin() and cos() in the loop entirely.

Well, it turns out that there is a way. The formula is presented in Numerical Recipes in C, pages 178-179, and takes advantage of the following recurrence relation:

cos (theta+delta) = cos(theta) - [alpha*cos(theta) + beta*sin(theta)]
sin (theta+delta) = sin(theta) - [alpha*sin(theta) - beta*cos(theta)]

Thus, given the previous values of cos(theta) and sin(theta), and the values alpha & beta, you can get the next values of sine and cosine with four multiplies and four adds. It is important to note that the values in the square brackets should be evaluated first in order to preserve significance of the result if delta is small.

The values alpha and beta are given as follows:

alpha = 2 * sin²(delta/2)
beta = sin(delta)

[In case your browser doesn't properly display ISO Latin 1, alpha is two times sine-squared of delta over two.]

Thus, the above loop would be recoded as follows:

double theta_increment = (2*PI) / nslices;
int    i = 0;
double Tcos = 1;   // Start from theta = zero.
double Tsin = 0;

double beta  = sin (theta_increment);
double alpha = sin (theta_increment/2);
alpha = 2 * alpha * alpha;

do
{   double Ncos, Nsin;   // New Values

    ...     // Use Tcos and Tsin

    Ncos = (alpha * Tcos) + (beta * Tsin);
    Nsin = (alpha * Tsin) - (beta * Tcos);

    Tcos -= Ncos;
    Tsin -= Nsin;

} while (++i < nslices);

Now here's the interesting part: this method preserves accuracy very well. In my test runs (90Mhz Pentium), I got the following results:

Increment (Degrees) Revolutions Max Error Sin() Max Error Cos() Standard Speed (sec) Incremental Speed (sec)
0.00001 1 3.496e-13 2.648e-13 83.8 9.36
0.001 100 1.820e-12 1.816e-12 87.3 9.29
0.1 10,000 4.113e-12 4.114e-12 90.3 9.21
1.0 100,000 6.828e-116.849e-11 90.49.50
10.0 1,000,000 8.1934e-10 8.1620e-10 94.29.34

That's phenominal -- that last case is the result of 36 million iterations, and is still accurate to nine decimal places! I was expecting a much faster degradation of the results, given that you continually accumulate rounding errors. In addition, I'd expected to be bitten by the finite precision of the sin() function to generate alpha and beta. Nevertheless, there is some small error in this method. It's easy enough to add a "resynch" counter to the loop to re-synchronize the Tcos and Tsin values (via the sin() and cos() functions) every N calculations if you feel that it's important to maintain accuracy indefinitely.


Steve Hollasch / Last Update 1998 December 16