A fast algorithm, called the fast Fourier transform on multipoles (FFTM) method, is developed for efficient solution of the integral equation in the boundary element method (BEM). This method employs the multipole and local expansions to approximate far field potentials, and uses the fast Fourier transform (FFT) to accelerate the multipole to local translation operator based on its convolution nature. The series of uncoupled convolutions allows further speed up in the algorithm through parallel computation. In this paper, we present the results of using the FFTM algorithm for solving large-scale three-dimensional electrostatic problems. It is demonstrated that the method can give accurate results with relatively low order of expansion. It is also found that the serial version of the algorithm has computational complexities of O(N/sup a/), where a ranges from 1.0 to 1.4 for computational time, and from 1.1 to 1.2 for memory storage requirement. Significant speedup is also observed in the parallel implementation of FFTM using up to 16 processors on an IBM-p690 supercomputer.