We describe a method for calculating the matrix elements of the Coulomb operator for Gaussian basis sets using an intermediate discrete Fourier transform of the density. Our goals are the same as those of the Gaussian and augmented-plane-wave method of Parrinello and co-workers [M. Krack and M. Parrinello, Phys. Chem. Chem. Phys. 2, 2105 (2000)], but our techniques are quite different. In particular, we aim at much higher numerical accuracy than typical programs using plane wave expansions. Our method is free of the effects of periodic images and yields full precision. Other low-scaling methods for the Coulomb operator are compared to the Fourier transform method with regard to numerical precision, asymptotic scaling with molecular size, asymptotic scaling with basis set size, onset point (the size of the calculation where the method outperforms traditional Gaussian integral techniques by a factor of 2), and the ability to calculate the Hartree-Fock exchange operator. The Fourier transform method is superior to alternatives by most criteria. In particular, for typical molecular applications it has an earlier onset point than fast multipole methods.