All formulas that are necessary for deriving not only upper (EU) but also lower (EL) variational bounds on the energy of systems featuring a few nonrelativistic particles are obtained with trial functions in the form of expansions in multidimensional Gaussian functions or exponentials. For potentials that are used most widely, all matrix elements are expressed in terms of known functions, a circumstance that simplifies considerably relevant numerical calculations. This is so for systems featuring an arbitrary number of particles in the case of a Gaussian basis and for three-particle systems in the case of an exponential basis. Numerical results for EU and EL, which are characterized by record accuracies, are presented for some Coulomb and nuclear systems such as the He atom; the e+e−e−, ppμ−, 3α, and 4α systems; and hypertritium (pnΛ). Lower bounds with exponential trial functions are obtained for the first time (the corresponding formulas are presented for the first time as well); for a Gaussian basis, lower bounds for Coulomb systems have not been known either. Given EL and EU, limits within which the exact value of energy, E0, lies can be indicated with confidence. Moreover, an analysis of the correlation between EL and EU with increasing number of terms in the expansion of the trial function makes it possible to improve the accuracy (at least by one order of magnitude) of the value E∞ extrapolated to infinity. By considering specific examples, it is shown that the exponential basis is advantageous in relation to the Gaussian one.