Using energy operators RC1-nRD1-m, RC1-nr12-m, and r12-nr13-m with small (n, m) values is fundamental in electronic structure calculations. Analytical integrations of the cases (n, m) = (1, 0) and (0, 1) are based on the Laplace transformation with the integrand exp(-a2t2), the other cases are based on the Laplace transformation with the integrand exp(-a2t) and the two-dimensional version of the Boys function. These analytic expressions, with Gaussian function integrands, are useful for manipulation with higher moments of interelectronic distances, for example, in correlation calculations. The equations derived help to evaluate the one-, two-, and three-electron Coulomb integrals, òρ(1)RC1-nRD1-mdr1, òρ(1)ρ(2)RC1-nr12-mdr1dr2, and òρ(1)ρ(2)ρ(3)r12-nr13-mdr1dr2dr3, wherein ρ(i) is the one-electron density describing the electron clouds in molecules, solids, or any media or ensemble of materials. Analytical solutions to integrals are more useful than numeric solutions; however, the former is not available in many cases. We evaluate these integrals numerically, even more so, the òf(ρ(1))dr1 to òf(ρ(1),ρ(2),ρ(3))dr1dr2dr3 with the analytical function f. For this task, the commonly used density functional theory numerical integration scheme has been elaborated to 6 and 9 dimensions via Descartes product. More importantly, this numerical integration scheme works not only for Gaussian type but also for Slaterian types. Analogy is commented on in terms of the powerful empirical correction between quantum potential energy correction and the empirically corrected Newton’s universal law of gravity in the explanation of dark matter and energy, as well as its relation to Hartree-Fock and Kohn-Sham formalisms.