Fast and highly accurate computation of Chebyshev expansion coefficients of analytic functions

Chebyshev expansion coefficients can be computed efficiently by using the FFT, and for smooth functions the resulting approximation is close to optimal, with computations that are numerically stable. Given sufficiently accurate function samples, the Chebyshev expansion coefficients can be computed to machine precision accuracy. However, the accuracy is only with respect to absolute error, and this implies that very small expansion coefficients typically have very large relative error. Upon differentiating a Chebyshev expansion, this relative error in the small coefficients is magnified and accuracy may be lost, especially after repeated differentiation. At first sight, this seems unavoidable. Yet, in this paper, we focus on an alternative computation of Chebyshev expansion coefficients using contour integrals in the complex plane. The main result is that the coefficients can be computed with machine precision relative error, rather than absolute error. This implies that even very small coefficients can be computed with full floating point accuracy, even when they are themselves much smaller than machine precision. As a result, no accuracy is lost after differentiating the expansion, and even the 100th derivative of an analytic function can be computed with near machine precision accuracy using standard floating point arithmetic. In some cases, the contour integrals can be evaluated using the FFT, making the approach both highly accurate and fast.