We present a public code to generate random fields with an arbitrary probability distribution function (PDF) and an arbitrary correlation function. The algorithm is cosmology-independent, applicable to any stationary stochastic process over a three dimensional grid. We implement it in the case of the matter density field, showing its benefits over the lognormal approximation, which is often used in cosmology for generation of mock catalogues. We find that the covariance of the power spectrum from the new fast realizations is more accurate than that from a lognormal model. As a proof of concept, we also apply the new simulation scheme to the divergence of the Lagrangian displacement field. We find that information from the correlation function and the PDF of the displacement-divergence provides modest improvement over other standard analytical techniques to describe the particle field in the simulation. This suggests that further progress in this direction should come from multi-scale or non-local properties of the initial matter distribution.