  Advances in photonics require photon-number resolved simulations of quantum
optical experiments with Gaussian states. We demonstrate a simple and versatile
method to simulate the photon statistics of general multimode Gaussian states.
The derived generating functions enable simulations of the photon number
distribution, cumulative probabilities, moments, and factorial moments of the
photon statistics of Gaussian states as well as of multimode photon-added and
photon-subtracted Gaussian states. Numerical results are obtained by automatic
differentiation of these generating functions by employing the software
framework PyTorch. Our approach is particularly well suited for practical
simulations of the photon statistics of quantum optical experiments in
realistic scenarios with low photon numbers, in which various sources of
imperfections have to be taken into account. As an example, we calculate the
detection probabilities for a recent multipartite time-bin coding quantum key
distribution setup and compare them with the corresponding experimental values.
