  In this paper, we provide an algorithm and general framework for the
simulation of photons passing through linear optical interferometers. Given $n$
photons at the input of an $m$-mode interferometer, our algorithm computes the
probabilities of all possible output states with time complexity
$O\left({n\binom{n+m-1}{m-1}}\right)$, linear in the number of output states
$\binom{n+m-1}{m-1}$. It outperforms the permanent-based method by an
exponential factor, and for the restricted problem of computing the probability
for one given output it improves the time complexity over the state-of-the-art
for the permanent of matrices with multiple rows or columns, with a tradeoff in
the memory usage. Our algorithm also has additional versatility by virtue of
its use of memorisation -- the storing of intermediate results -- which is
advantageous in situations where several input states may be of interest.
Additionally it allows for hybrid simulations, in which outputs are sampled
from output states whose probability exceeds a given threshold, or from a
restricted set of states. We consider a concrete, optimised implementation, and
we benchmark the efficiency of our approach compared to existing tools.
