When comparing approximate Gaussian process (GP) models, it can be helpful to be able to generate data from any GP. If we are interested in how approximate methods perform at scale, we may wish to generate very large synthetic datasets to evaluate them. Na\"{i}vely doing so would cost (\order{n^3}) flops and (\order{n^2}) memory to generate a size (n) sample. We demonstrate how to scale such data generation to large (n) whilst still providing guarantees that, with high probability, the sample is indistinguishable from a sample from the desired GP.