Multiplication of polynomials in Chebyshev basis For polynomials in the monomial basis like $p_n(x) = \sum_{k=0}^N a_k x^k $, the product of 2 polynomials is can be either found though the convolution of the 2 corresponding polynomial vectors or with FFT/IFFT.
I wonder, if there exists a "numerical recipe" to compute the product of 2 polynomials like $ p_n(x) = \sum_{k=0}^N c_k T_k(x)$ (i.e. represented in the Chebyshev basis).
 A: Yes, they boil down to sums of things that look just like convolution or correlation of the coefficients. See equations 2.6, 2.7, and 2.8 here for the exact formulas:
http://www2.mathematik.hu-berlin.de/~gaggle/S09/AUTODIFF/projects/papers/baszenski_fast_polynomial_multiplication_and_convolutions_related_to_the_discrete_cosine_transform.pdf
Furthermore, the sums can be efficiently calculated using FFT-related techniques, which also seems to be outlined in that paper, though I didn't look at the details.
Edit: To make this a complete answer that does not require visiting the reference, here are the relevant formulas:
For $n\in\mathbb{N}$, $p_n$ and $q_n$ are $n^\text{th}$ order polynomials written as Chebyshev series in the $T_k$ polynomials:
$$
p_n = \frac{a_0}{2}+\sum_{k=1}^n a_k T_k,\;\;\;q_n = \frac{b_0}{2}+\sum_{k=1}^n b_k T_k
$$
The coefficients $a_k$ and $b_k$ are real. Then the product of the two polynomials has order $2n$ and is denoted $r_{2n}$. The Chebyshev series form of that product is
$$
r_{2n} = \frac{c_0}{2}+\sum_{k=1}^{2n} c_k T_k
$$
And the $c_k$ coefficients can be calculated from the original $a_k$ and $b_k$ coefficients as:
$$
2c_k = \left\{
\begin{array}{lcl}
\displaystyle a_0b_0 + 2\sum_{l=1}^n a_l b_l, & {} & k=0 \\
\displaystyle \sum_{l=0}^k a_{k-l}b_l + \sum_{l=1}^{n-k}\left(a_l b_{l+k}+a_{l+k}b_l\right), & {} & k=1,\ldots,n-1 \\
\displaystyle \sum_{l=n-k}^n a_{k-l} b_l, & {} & k=n,\ldots,2n
\end{array}
\right.
$$
