Fast evaluation of sums of products of two sinusoidal functions

Question

$s_1, s_2,\ldots,s_N$ is a set of real numbers. It is required to evaluate the sums \begin{eqnarray} S_1(n_1,n_2)&=&\sum_{j=1}^N \cos(n_1 s_j)\cos(n_2 s_j)\nonumber\\ S_2(n_1,n_2)&=&\sum_{j=1}^N \cos(n_1 s_j)\sin(n_2 s_j)\nonumber\\ S_3(n_1,n_2)&=&\sum_{j=1}^N \sin(n_1 s_j)\sin(n_2 s_j)\nonumber \end{eqnarray} for $n_1, n_2=1,2,3\ldots,M$. The length $N$ of ${\bf s}$ is typically of the order of several hundred to a thousand while $M$ is several tens of thousands. A fast numerical algorithm is required. For small $M$ the three matrices $S_k(n_1,n_2)$ can quickly be calculated by matrix multiplication, but this is not a realistic option (at least not for my PC) when M approaches $10^5$. In my application it is not necessary to store the $S_k(n_1,n_2)$ -- these can used and then discarded as they are calculated -- but the values are needed twice, at different stages of a larger algorithm, hence an efficient calculation algorithm is important.