Summary Currently, gravimetric forward modelling of mass density structures with arbitrary geometries and density distributions typically involves subdividing the mass body into individual geometric elements (such as rectangular prisms), calculating their gravitational contributions that are then summed up to obtain the gravitational attraction of the whole body. To achieve a more accurate approximation of the true geometric shape and density distribution, this rectangular prism model requires fine dividing, which significantly increases computational load and reduces numerical efficiency. To address this issue, we propose the algorithm for gravimetric forward modeling of arbitrary geometric shapes and density distributions in spectral domain that significantly improves numerical efficiency while preserves computational accuracy. The novelty of our proposed algorithm lies in dividing the masses into multiple layers of equal thickness in the vertical direction, providing constant upper and lower bounds. This allows to extended Parker's formulas and apply the Fast Fourier Transform (FFT) to increase numerical efficiency. The algorithm is tested using synthetic models and then used to compute gravitational effects of topography and sediments using real data from Tibet. Results show high accuracy and numerical efficiency than rectangular prism approach.