This paper presents a method for the computation of polar harmonic transforms that is fast and efficient. The method is based on the inherent recurrence relations among harmonic functions that are used in the definitions of the radial and angular kernels of the transforms. The employment of these relations leads to recursive strategies for fast computation of harmonic function-based kernels. Polar harmonic transforms were recently proposed and have shown nice properties for image representation and pattern recognition. The proposed method is 10-time faster than direct computation and five-time faster than fast computation of Zernike moments.