Fast and accurate numerical modeling of gravity and magnetic anomalies is the basis of field-data inversion and quantitative interpretation. In gravity and magnetic prospecting, the computation and memory requirements of practical modeling is still a significant issue, which leads to the difficulty of using efficient and detailed inversions for large-scale complex models. A new 3D numerical modeling method for gravity and magnetic anomaly in a mixed space-wavenumber domain is proposed to mitigate the difficulties. By performing a 2D Fourier transform along two horizontal directions, 3D partial differential equations governing gravity and magnetic potentials in the spatial domain are transformed into a group of independent 1D differential equations wrapped with different wavenumbers. Importantly, the computation and memory requirements of modeling are greatly reduced by this method. A modeling example with 4,040,100 observations can be finished in approximately 28 s on a desktop using a single core, and the independent differential equations are highly parallel among different wavenumbers. The method preserves the vertical component in the space domain, and thus a mesh for modeling can be finer at a shallower depth and coarser at a deeper depth. In general, the new method takes into account the calculation accuracy and the efficiency. The finite-element algorithm combined with a chasing method is used to solve the transformed differential equations with different wavenumbers. In a synthetic test, a model with prism-shaped anomalies is used to verify the accuracy and efficiency of the proposed algorithm by comparing the analytical solution, our numerical solution, and a well-known numerical solution. Furthermore, we have studied the balance between computational accuracy and efficiency using a standard fast Fourier transform (FFT) method with grid expansion and the Gauss-FFT method. A model with topography is also used to explore the ability of modeling topography with our method. The results indicate that the proposed method using the Gauss-FFT method has characteristics of fast calculation speed and high accuracy.