Computational Science Stack Exchange is a question and answer site for scientists using computers to solve scientific problems. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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
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

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

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

Your Answer


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.