  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.
