Abstract We propose a new modified Parker's method for efficient gravitational forward modeling and inversion using general polyhedral models. We have made several important modifications to the classical method, including: (a) The new method is now applicable to a general polyhedron represented by triangulated surface or tetrahedral mesh, and with arbitrarily variable 3D density distribution. (b) An optimal Fourier‐domain sampling strategy is used to improve the numerical accuracy of the new algorithm significantly. (c) A simple and effective automatic layering technique is introduced to accelerate the convergence rate of Parker's method. The method is demonstrated using both synthetic and real polyhedral models, including a sphere model approximated by a polyhedron, two asteroids, a digital elevation model in the Himalaya region, and the Yucca Flat basin model in Nevada. The numerical results show that, compared with analytical solutions of polyhedral models in the space domain, the modified Parker's method can improve the computational efficiency by several orders of magnitude while obtaining almost the same simulation results. The difference is well below existing instrumentation error level. By embedding the new forward algorithm into an iterative process, it can be used for fast inversion of density interfaces. Our new method is suitable for the efficient modeling and inversion of gravitational potential, gravitational vector, and gravitational gradient tensor caused by polyhedral models with a large number of faces, representing geological abnormal bodies, asteroids, and single or multilayer density interface models with triangulated surfaces.