I stumbled upon this question and I think I have an interesting answer. I will make use of the FFT, please remark that in this argument I do not deal with the numerical error of the FFT but there are explicit error bounds on the FFT so they may be taken into account (you can find them in Higham, Accuracy and Stability...); moreover I will use the representation of Fourier Series with exponentials because it makes the argument easier.
I suppose we are on $[0,1]$.
To compute the FFT we evaluate the trigonometric polynomial at $N$ equispaced points and run the Cooley-Tuckey algorithm. The algorithm returns us $X_k$ such that $$ X_k = \sum_{p=-\infty}^{+\infty} \hat{f}(k+pN) $$ where $\hat{f}(k)$ is the $k$-th Fourier coefficient of $f$ (with respect to the exponential basis).
Since $f$ is a trigonometric polynomial $\hat{f}(k)=0$ for $k>K$, so there exists an $N$ such that $X_k = \hat{f}(k)$.
We observe now that $\hat{f'}(k)=2\pi i k \hat{f}(k)$. Therefore,
$$||f'||_{\infty}\leq 2\pi K \sum_{-K}^K |X_k|. $$
We now run the inverse FFT of size M on the $X_k$. This gives us the value of the trigonometric polynomial at $M$ equispaced points, call them $x_1, \ldots, x_M$. Then $$ \max_{i=1,\ldots, M} f(x_i)\leq \max f(x)\leq \max_{i=1,\ldots, M} f(x_i)+ 2\pi \frac{K}{M} \sum_{-K}^K |X_k|. $$
The FFT and Inverse FFT are really fast, they are computed in O(N log(N)). Moreover, if you already know the coefficients of the trigonometric polynomial in the $\sin$, $\cos$ basis you can convert them easily to the exponential basis.