The author has previously extended the theory of regular and irregular primes to the setting of arbitrary totally real number fields. It has been conjectured that the Bernoulli numbers, or alternatively the values of the Riemann zeta function at odd negative integers, are evenly distributed modulo p for every p. This is the basis of a well-known heuristic, given by Siegel, for estimating the frequency of irregular primes. So far, analyses have shown that if Q(\sqrt{D}) is a real quadratic field, then the values of the zeta function \zeta_{D}(1-2m)=\zeta_{Q(\sqrt{D})}(1-2m) at negative odd integers are also distributed as expected modulo p for any p. However, it has proven to be very computationally intensive to calculate these numbers for large values of m. In this paper, we present the alternative of computing \zeta_{D}(1-2m) for a fixed value of D and a large number of different m.