Connections between recurrence relations, spectral theory, and the approximation theory of Laurent series lead to efficient, numerically stable algorithms for the analysis and synthesis of linear combinations of special functions. These include algorithms both for computing the coefficients in linear combinations of the functions, given the values of these linear combinations at certain points, and, vice versa, for evaluating such linear combinations at those points, given the coefficients in the linear combinations. The costs of the algorithms for a linear combination of n functions are bounded by a constant times n

log(n) at any fixed precision; moreover, the constant is reasonably small in single- or double-precision arithmetic, and is the same for all families of special functions satisfying three-term recurrence relations, including the classical orthogonal polynomials. Inter alia, this yields fast spherical harmonics transforms, as used for large-scale weather and climate simulations across the globe.