  In a recent work we presented a recursive algorithm to compute the matrix
elements of a generic Gaussian transformation in the photon-number basis. Its
purpose was to evolve a quantum state by building the transformation matrix and
subsequently computing the matrix-vector product. Here we present a faster
algorithm that computes the final state without having to generate the full
transformation matrix first. With this algorithm we bring the time complexity
of computing the Gaussian evolution of an $N$-dimensional $M$-mode state from
$O(MN^{2M})$ to $O(M(N^2/2)^M)$, which is an exponential improvement in the
number of modes. In the special case of high squeezing, the evolved state can
be approximated with complexity $O(MN^{M})$. Our new algorithm is
differentiable, which means we can use it in conjunction with gradient-based
optimizers for circuit optimization tasks. We benchmark our algorithm by
optimizing circuits to produce single photons, Gottesman-Kitaev-Preskill states
and NOON states, showing that it is up to one order of magnitude faster than
the state of the art.
