  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.
