Abstract. The spectral method with discrete spherical harmonics transform plays an important role in many applications. In spite of its advantages, the spherical harmonics transform has a drawback of high computational complexity, which is determined by that of the associated Legendre transform, and the direct computation requires time of O(N 3 ) for cut-off frequency N . In this paper, we propose a fast approximate algorithm for the associated Legendre transform. Our algorithm evaluates the transform by means of polynomial interpolation accelerated by the Fast Multipole Method (FMM). The divideand-conquer approach with split Legendre functions gives computational complexity O(N 2 log N ). Experimental results show that our algorithm is stable and is faster than the direct computation for N ≥ 511.