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-11 | 6.849e-11 | 90.4 | 9.50 |

10.0 | 1,000,000 | 8.1934e-10 | 8.1620e-10 | 94.2 | 9.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