Thursday, May 23, 2024

[ftsgsdrj] rapid evaluation of sine and cosine via recurrence

evenly spaced evaluations of sine and cosine are common, for example, in Fourier transforms or plotting circles.  the following recurrence is given in the first few paragraphs section 5.5 "Recurrence Relations and Clenshaw's Recurrence Formula" of "Numerical Recipes" Second Edition (section 5.4 in the Third Edition, because they have switched to indexing book sections (and C arrays) from zero).

let
alpha = 2*(sin(delta/2))^2
beta = sin(delta)

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

evaluate the expression in square brackets first to avoid losing significance; do NOT use associativity to move the square brackets.

I think the second formula is equivalent via commutativity to the following, but I'm not completely sure with respect to speed and numerical stability.  writing it the following way also makes it more obvious that the expressions in square brackets are matrix-vector multiplication with a constant matrix [alpha , beta ; -beta , alpha], reminiscent of a rotation matrix.  maybe you have fast 2x2 matrix-vector multiplication.

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

Because it is a recurrence, use the new values of cos(theta + delta) and sin(theta + delta) to calculate the next increment cos(theta + delta + delta) and sin(theta + delta + delta ), and so forth.

future work: evaluate the performance and accuracy of the recurrence.  does it work best when delta is small?

this post was motivated by not being able to find the recurrence in Wikipedia, nor other neat or even standard tricks for computing sine and cosine.  the above recurrence is not given on https://en.wikipedia.org/w/index.php?title=Sine_and_cosine&oldid=1204203183#Software_implementations nor https://en.wikipedia.org/w/index.php?title=Trigonometric_functions&oldid=1204409434 , though the latter article seems to currently be being aggressively edited and regularly trimmed by a highly opinionated editor, so who knows what useful information has been deleted.

the other trick that I know of is repeatedly halving the angle to make the Maclaurin (Taylor) series converge faster, then using double-angle formulae to repeatedly double the angle back to the original angle.

No comments :