The density inversion of gravity data is commonly achieved by discretizing the subsurface into prismatic cells and calculating the density of each cell. During this process, a weighting function is introduced to the iterative computation to reduce the skin effect during the inversion. Thus, the computation process requires a significant number of matrix operations, which results in low computational efficiency. We have adopted a density inversion method with nonlinear polynomial fitting (NPF) that uses a polynomial to represent the density variation of prismatic cells in a certain space. The computation of each cell is substituted by the computation of the nonlinear polynomial coefficients. Consequently, the efficiency of the inversion is significantly improved because the number of nonlinear polynomial coefficients is less than the number of cells used. Moreover, because representing the density change of all of the cells poses a significant challenge when the cell number is large, we adopt the use of a polynomial to represent the density change of a subregion with fewer cells and multiple nonlinear polynomials to represent the density changes of all prism cells. Using theoretical model tests, we determine that the NPF method more efficiently recovers the density distribution of gravity data compared with conventional density inversion methods. In addition, the density variation of a subregion with 8 × 8 × 8 prismatic cells can be accurately and efficiently obtained using our cubic NPF method, which also can be used for noisy data. Finally, the NPF method is applied to real gravity data in an iron mining area in Shandong Province, China. Convergent results of a 3D perspective view and the distribution of the iron ore bodies are acquired using this method, demonstrating the real-life applicability of this method.