Large scale quantum chemistry calculations often involve evaluating and processing millions of Gaussian integrals. Whereas literature in this field,??like other areas of scientific computing, has focused on devising more efficient algorithms for evaluating these integrals, relatively little has been done to fathom the effects numerical errors may have on their accuracy. In this work we present several methods for computing rigorous bounds on the incomplete gamma function, known in this context as $F_m(T)$, which is a core quantity used in Gaussian integral evaluation. Based on a computational paradigm called interval arithmetic, these bounds are guaranteed to contain the true numerical result of $F_m(T)$, notwithstanding the presence of rounding and truncation errors. Experimental results are analysed to determine the best numerical approaches for ensuring that the bounds are not only rigorous, but also sufficiently precise. A worst-case error analysis of existing numerical techniques is also presented as a consequence of this work.