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 to do 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.


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.

  • $\begingroup$ 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? $\endgroup$
    – Victor Liu
    Nov 28 '12 at 3:20
  • $\begingroup$ 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. $\endgroup$
    – user335
    Nov 28 '12 at 3:22
  • $\begingroup$ 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. $\endgroup$
    – Erik P.
    Nov 28 '12 at 19:15
  • $\begingroup$ You're right. My answer was based on the assumption that all intermediate values ($1,...,m-1$) will also be used. $\endgroup$
    – user335
    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)$.


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 exceptionally well, especially in the case where $\varepsilon$ is tiny.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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