Optimization of harmonics calculation I need to compute the sine and cosine of an argument along with n &quot;harmonics&quot;
\begin{matrix} \sin(x) &amp; \cos(x) \\
  \sin(2x) &amp; \cos(2x) \\
\cdots \\
\sin(nx) &amp; \cos(nx)
\end{matrix}
This takes a lot of computing time and since this is a real-time system, I need to optimize this calculation.
Here's my solution so far :
I defined a rotation matrix $$
R = \begin{pmatrix}
    \cos(x) \ -\sin(x) \\
    \sin(x) \ \cos(x) \\
    \end{pmatrix}
$$
I compute $[\cos(nx); \sin(nx)] = R  [\cos((n-1)x); \sin((n-1)x)]$; I simply iterate to compute the missing $n-1$ harmonics.
This method is 3 times fast than computing all harmonics separately.
I'm curious, are there better solutions? I know a look-up table could be even faster but I would prefer not to go this way.
No, you can't do it much simpler or faster than that.  In the old days, subscripts used to cost time, not sure that is true any more.

#include &lt;math.h&gt;
#include &lt;stdio.h&gt;

//=========================================================
int main( int argCount, char *argValues[] )
{
        double x = 1.2;

        double C = cos( x );
        double S = sin( x );
        
        double CA[10];
        double SA[10];
        
        CA[0] = 1.0;
        SA[0] = 0.0;

        CA[1] = C;
        SA[1] = S;
        
        for( int h = 2; h &lt; 10; h++ )
        {
            CA[h] = CA[h-1] * C - SA[h-1] * S;
            SA[h] = CA[h-1] * S + SA[h-1] * C;
            
            printf( "%2d   %10.6f %10.6f      %10.6f %10.6f\n",
                    h, CA[h], SA[h], cos( h * x ), sin( h * x ) );
        }
        
        return 0;
}
//=========================================================

Results:

 2    -0.737394   0.675463       -0.737394   0.675463
 3    -0.896758  -0.442520       -0.896758  -0.442520
 4     0.087499  -0.996165        0.087499  -0.996165
 5     0.960170  -0.279415        0.960170  -0.279415
 6     0.608351   0.793668        0.608351   0.793668
 7    -0.519289   0.854599       -0.519289   0.854599
 8    -0.984688  -0.174327       -0.984688  -0.174327
 9    -0.194330  -0.980936       -0.194330  -0.980936

This implementation is based on complex multiplication rather than matrix rotation.  In this case, conceptually they are different, mechanically they are identical.
You will accumulate a little bit of error, but for small sets this remains in the $10^{-17}$ region.
See:

Angle Addition Formulas from Euler's Formula