  We introduce an exact classical algorithm for simulating Gaussian Boson
Sampling (GBS). The complexity of the algorithm is exponential in the number of
photons detected, which is itself a random variable. For a fixed number of
modes, the complexity is in fact equivalent to that of calculating output
probabilities, up to constant prefactors. The simulation algorithm can be
extended to other models such as GBS with threshold detectors, GBS with
displacements, and sampling linear combinations of Gaussian states. In the
specific case of encoding non-negative matrices into a GBS device, our method
leads to an approximate sampling algorithm with polynomial runtime. We
implement the algorithm, making the code publicly available as part of Xanadu's
The Walrus library, and benchmark its performance on GBS with random Haar
interferometers and with encoded Erd\H{o}s-Renyi graphs.
