  The statistical translation invariance of cosmological random fields is
broken by a finite survey boundary, correlating the observable Fourier modes.
Standard methods for generating Gaussian fields either neglect these
correlations, or are costly, or both. Here we report on a fast and exact
simulation method applicable to a wide class of two-point statistics that
requires the simulation of a periodic grid of only twice the survey side with
fast Fourier transforms. Super-survey modes, dominating the covariance of power
spectra beyond linear scales in galaxy surveys and causing the correlation of
large and small scales, "beat coupling", or "super-sample" covariance, are
precisely accounted for in non-linear transformations of the Gaussian field. As
an application, we simulate the CFHTLS $\sim 7^\circ \times 7^\circ$ W1 galaxy
density field, modeled as a Poisson sampling of a lognormal density field. We
show that our simulations produce power spectra, $A^*$-power spectra,
counts-in-cells probability distributions as well as covariances perfectly
consistent with the data. In addition, our technique reproduces the information
plateau beyond linear scales as observed previously in Sloan Digital Sky Survey
galaxy catalogs and in $N$-body simulations. Our method is thus an efficient
yet powerful simulation and prediction tool for galaxy survey data and
covariances.
