We introduce a new method to price American options based on Chebyshev interpolation. In each step of a dynamic programming time-stepping we approximate the value function with Chebyshev polynomials. The key advantage of this approach is that it allows us to shift the model-dependent computations into an offline phase prior to the time-stepping. In the offline part a family of generalized (conditional) moments is computed by an appropriate numerical technique such as a Monte Carlo, PDE, or Fourier transform based method. Thanks to this methodological flexibility the approach applies to a large variety of models. Online, the backward induction is solved on a discrete Chebyshev grid, and no (conditional) expectations need to be computed. For each time step the method delivers a closed form approximation of the price function along with the options' delta and gamma. Moreover, the same family of (conditional) moments yield multiple outputs including the option prices for different strikes, maturities, and different payoff profiles. We provide a theoretical error analysis and find conditions that imply explicit error bounds for a variety of stock price models. Numerical experiments confirm the fast convergence of prices and sensitivities. An empirical investigation of accuracy and runtime also shows an efficiency gain compared with the least-squares Monte Carlo method introduced by Longstaff and Schwartz [Rev. Financ. Stud., 14 (2001), pp. 113-147]. Moreover, we show that the proposed algorithm is flexible enough to price barrier and multivariate barrier options.