We present a fast algorithm based on polynomial interpolation to approximate the coefficient matrix arising from multi-term fractional diffusion equations, which is discretized by the implicit finite difference scheme with the shifted Gr¨unwald formula. The approximate matrix can be constructed in O(N) operations and requires O(N) storage, and costs only O(N log N) complexity for the matrix-vector multiplication, where N is the number of grid points. A multigrid method is proposed to solve the approximation system. Finally we extend our algorithm to solve distributed order space fractional differential equations, which is approximated by a multi-term fractional differential equation. Numerical results are given to demonstrate the accuracy and efficiency of the proposed algorithm.