How should I go about calculating sines or cosines without the errors introduced solely by the transcendence of pi?

Question

This may be irrelevant for many people but I'm asking this question because I don't need to do these things in real time, I don't care if modifying a 1-minute sound recording takes 5 minutes to process. My question has two parts.

#1. Conventional implementations of sine or cosine functions expect the argument to be in radians, which has one serious drawback. Many people will tell me to first reduce the range of the argument. But $\pi$ is transcendental, right? Therefore, if I try to represent $5 \pi$ in IEEE doubles as accurately as I can and if I then try to represent $\pi$ in IEEE doubles as accurately as I can, dividing the two approximations will not give me $5$. So I can hardly subtract $2 \pi$ twice without introducing further errors by doing this. Also, if I wanted to output the sine or cosine values themselves as IEEE doubles, I'd like the algorithm to be able to calculate them in such a way that the error would be unmeasurable with IEEE doubles. But this requirement can never be met if the approximation of $\frac{\pi}{3}$ divided by the approximation of $\frac{\pi}{9}$ differs from $3$ so much that the error is noticeable in IEEE doubles.

#2. Let's say that I only want to find sines or cosines of some rational parts of $\pi$. For this, it might be meaningful to use Chebyshev polynomials and calculate some insanely precise lookup tables (maybe some double-doubles or whatever) before launching any actual signal processing code. But let's say I want to find $2 \sin \left( \frac{\pi}{9} \right)$ and the corresponding equation goes like this:

$$ 3 x - x^{3} = \sqrt{3} \tag{1} \label{1}$$

Even if I knew about an excellent algorithm for finding cube roots, I could not use it here because that algorithm would most probably be meant for finding cube roots of real numbers, not of complex numbers. And if I were supposed to convert the complex numbers to polar form and back again, I would be introducing even more errors into this. So my idea would be to somehow find roots of Chebyshev polynomials without involving the conversion to polar form. But so far I only know about one person who seems to have described a similar approach in great detail and I'm not sure if anyone else has (if you're interested, I can post the link to the paper).

So my question is: What approach[es] would you suggest me to take if I only care about accuracy and don't care about the time taken? I mean, as I've said, I don't need to store the final data with an accuracy of some $100$ decimal digits, but IEEE doubles would be good.

Answer 1

Most things that you do in a computer have finite precision. The approach to dealing with this is straight forward:

1. Determine the precision that your application needs
2. Analyze error propagation and accumulation in your algorithm
3. Choose a data type that has enough precision to accommodate your target precision PLUS the numerical error.
4. Run algorithm using your chosen data type(s) precision and convert to your target data-type/precision.

Simple example: For most audio applications, 120 dB SNR would be considered as "good enough". If you have an algorithm than adds about 10 dB of numerical noise (through rounding etc), you want to use a data type with at least 130 dB SNR. IEEE single precision floating point has 144 dB SNR and is hence a very popular choice for audio signal procession.

In your specific example: IEEE double precision gives you about 324 dB of SNR. If that's not good enough than use a higher precision data type and convert back when you done. Double Double and Quad Double libraries are readily available (see for example https://github.com/scibuilder/QD)