Take the 2-minute tour ×
Computational Science Stack Exchange is a question and answer site for scientists using computers to solve scientific problems. It's 100% free, no registration required.

When evaluating cylindrical harmonics, one needs to evaluate trigonometric functions $\cos(m\theta)$ and $\sin(m\theta)$, potentially for large integer $m$ and $\theta\in[-\pi,\pi]$. What is the best way of doing this in C code? Currently, I just evaluate at the angle $m\theta$ but I would suspect that standard libraries lose accuracy at large arguments. I was considering using double angle formulas and the like to recursively reduce the magnitude of the arguments, but I'm wondering if that ends up incurring more error.

share|improve this question
add comment

3 Answers 3

up vote 3 down vote accepted

If you do it iteratively, by computing $\sin(n \theta)$ and $\cos(n \theta)$ from $\sin((n - 1) \theta)$ and $\cos((n - 1) \theta)$, then floating point errors will not blow up by accumulation. This is because the transition matrix is orthonormal. It is a simple rotation matrix.

share|improve this answer
    
Hm, that's a good point. However, even if the error's don't blow up, do I have some assurance that the result hasn't drifted away from what the "exact" answer should be? –  Victor Liu Nov 28 '12 at 3:20
    
I doubt that the results will be very different compared to the exact answer. You can try for $m = 1000$ and $\theta = \pi/m$ so you know what to expect. –  Chetan Jhurani Nov 28 '12 at 3:22
    
Whenever $n$ is even, I would use the double angle formula to cut the exponent in half. This corresponds to computing the $n$th power of the rotation matrix with the $O(\log(n))$ algorithm. –  Erik P. Nov 28 '12 at 19:15
    
You're right. My answer was based on the assumption that all intermediate values ($1,...,m-1$) will also be used. –  Chetan Jhurani Nov 28 '12 at 23:49
add comment

I've encountered a similar problem before, but I was more worried about speed than accuracy (see paper here). If your angle $\vartheta$ is the result of an $\arccos(\cdot)$, which is often the case in geometric computations, you can use Chebyshev polynomials which are defined as

$$T_k(x) = \cos(k\arccos(x)), \quad \mbox{or} \quad T_k(\cos(\vartheta)) = \cos(k\vartheta)$$

and can be evaluated quickly ans stably using the three-term recurrence relation

$$ T_k(x) = 2xT_{k-1}(x) - T_{k-2}(x), \quad T_0(x) = 1, \quad T_1(x) = x. $$

Thus, for your problem, you can evaluate $\cos(m\vartheta)$ in $m+1$ multiplications and additions, provided you have $\cos(\vartheta)$.

share|improve this answer
add comment

The procedure I prefer, which I have also seen used by some FFT implementations, also involves a three-term recurrence, but is a bit more stable than the computation of $\cos\,kx$ and $\sin\,kx$ from previous predecessors. (There is a nice analysis of these recurrences in Bulirsch/Stoer; see examples 2 and 3. A geometric way of looking at the instability is that although you trace out an approximate circle at first, in the long run you literally spiral out of control.)

The recurrence (effectively the same as example 4 in Bulirsch/Stoer, see that for a detailed analysis) proceeds like so:

$$\begin{align*} \cos(\theta+\varepsilon)&=\cos\,\theta-(p\cos\,\theta+q\sin\,\theta)\\ \sin(\theta+\varepsilon)&=\sin\,\theta-(p\sin\,\theta-q\cos\,\theta) \end{align*}$$

where we have the precomputed constants

$$p=2\sin^2\frac{\varepsilon}{2},\qquad q=\sin\,\varepsilon$$

As noted, they perform excpetionally well, especially in the case where $\varepsilon$ is tiny.

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.