We use first-order perturbation theory to represent the three-dimensional (3-D) seismic structure of the Earth. The background structure is assumed one-dimensional (1-D) with variations in only vertical direction. The perturbations are assumed 3-D volumetric inclusions with certain velocity contrast to the embedding structure. A convolutional type integral equation employing the Green's function of the 1-D medium is utilized to solve the total wavefield at any point in the 3-D medium. The 3-D perturbations acting as secondary sources are allowed to have arbitrary shapes and sizes, and the integral equation is solved using the numerical techniques. A flat Earth structure is assumed in the mathematical formulation and the spherical-to-flat transforms are employed to account for the Earth's sphericity. The 1-D Green's function is extracted using the stable reflectivity method. The integral equation consists of five integrals in the spectral and spatial domains. The Multilevel Fast Multipole Method (MLFMM) employed for the two integrals on the horizontal plane greatly helps reduce the computational time. The speed gain is characterized as logarithmic versus exponential if a traditional technique instead of the MLFMM is applied. The difference between the logarithmic and exponential time complexities becomes evident especially when the data size enlarges.