  We propose efficient classical algorithms which (strongly) simulate the
action of bosonic linear optics circuits applied to superpositions of Gaussian
states. Our approach relies on an augmented covariance matrix formalism to keep
track of relative phases between individual terms in a linear combination. This
yields an exact simulation algorithm whose runtime is polynomial in the number
of modes and the size of the circuit, and quadratic in the number of terms in
the superposition. We also present a faster approximate randomized algorithm
whose runtime is linear in this number. Our main building blocks are a formula
for the triple overlap of three Gaussian states and a fast algorithm for
estimating the norm of a superposition of Gaussian states up to a
multiplicative error. Our construction borrows from earlier work on simulating
quantum circuits in finite-dimensional settings, including, in particular,
fermionic linear optics with non-Gaussian initial states and Clifford
computations with non-stabilizer initial states. It provides algorithmic access
to a practically relevant family of non-Gaussian bosonic circuits.
