  When sampling the distribution P(phi) ~ exp(-|A phi|^2), a global heatbath
normally proceeds by solving the linear system A phi = eta, where eta is a
normal Gaussian vector, exactly. This paper shows how to preserve the
distribution P(phi) while solving the linear system with arbitrarily low
accuracy. Generalizations are presented.
