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 runI will work with the Cooley-Tuckey algorithmcomplex exponential basis. The algorithm returns us $X_k$ such that $$ X_k = \sum_{p=-\infty}^{+\infty} \hat{f}(k+pN) $$ where Let $\hat{f}(k)$ isbe the $k$-th Fourier coefficientcoefficients of $f$ (withthe Fourier series with respect to the exponential basis)basi.
Since Since $f$ is a trigonometric polynomial $\hat{f}(k)=0$ for $k>K$, so there exists ana $N$$K$ such that $X_k = \hat{f}(k)$ $\hat{f}(k)=0$ for $k>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|. $$$$||f'||_{\infty}\leq 2\pi K \sum_{-K}^K |\hat{f}(k)|. $$
We now run the inverse FFT of size M bigger than K on the $X_k$$\hat{f}(k)$. This The inverse Fast Fourier Transform using the Cooley-Tuckey algorithm is fast and numerically well behaved and gives us the value of the trigonometric polynomial at $M$ equispaced points, call them $x_1, \ldots, x_M$. Then
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|. $$$$ \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 |\hat{f}(k)|. $$
The FFT and Inverse FFT are really fast, theyThere are computed in O(N log(N)). Moreover, if you already know the coefficients of the trigonometric polynomialsome normalizations involved in the $\sin$FFT, $\cos$ basis you can convert them easily tobut the exponential basisargument should work.