  We develop closed-form expressions for the moments and cumulants of Gaussian
states when measured in the photon-number basis. We express the photon-number
moments of a Gaussian state in terms of the loop Hafnian, a function that when
applied to a $(0,1)$-matrix representing the adjacency of a graph, counts the
number of its perfect matchings. Similarly, we express the photon-number
cumulants in terms of the Montrealer, a newly introduced matrix function that
when applied to a $(0,1)$-matrix counts the number of Hamiltonian cycles of
that graph. Based on these graph-theoretic connections, we show that the
calculation of photon-number moments and cumulants are $#P$-hard. Moreover, we
provide an exponential time algorithm to calculate Montrealers (and thus
cumulants), matching well-known results for Hafnians. We then demonstrate that
when a uniformly lossy interferometer is fed in every input with identical
single-mode Gaussian states with zero displacement, all the odd-order cumulants
but the first one are zero. Finally, we employ the expressions we derive to
study the distribution of cumulants up to the fourth order for different input
states in a Gaussian boson sampling setup where $K$ identical states are fed
into an $\ell$-mode interferometer. We analyze the dependence of the cumulants
as a function of the type of input state, squeezed, lossy squeezed, squashed,
or thermal, and as a function of the number of non-vacuum inputs. We find that
thermal states perform much worse than other classical states, such as squashed
states, at mimicking the photon-number cumulants of lossy or lossless squeezed
states.
