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.