Optimization of harmonics calculation
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