Prismatic bodies with variable density contrasts are very useful in modeling and inversion of gravity anomalies for sedimentary basins, in which some simple mathematical functions offer better approximations for the density-depth relationship of the basin fill than a constant density model. We have explored Fourier-domain modeling of vector and tensor gravity anomalies due to 2D and 3D prismatic bodies following a wide class of variable density functions. We placed the emphasis on general polynomial models because they provided a flexible way to approximate arbitrary density functions. The recently proposed Gauss-fast Fourier transform method was applied to improve the accuracy of inverse Fourier transform. The numerical performance of the presented method was examined by comparing with space-domain analytical, seminumerical, or completely numerical solutions. Singularity and numerical stability of the algorithm were also analyzed with ranges of numerical stability for polynomial density models of different orders specified. Truncation errors of Fourier forward methods were estimated quantitatively using relevant model parameters, and the behavior was observed through theoretical analysis and model tests. Synthetic models, of 2D and 3D prisms with variable density contrasts, found the accuracy, efficiency, and flexibility of the method. We also interpreted a real data example of the gravity profile across the San Jacinto graben.