  Evaluating the Torontonian function is a central computational challenge in
the simulation of Gaussian Boson Sampling (GBS) with threshold detection. In
this work, we propose a recursive algorithm providing a polynomial speedup in
the exact calculation of the Torontonian compared to state-of-the-art
algorithms. According to our numerical analysis the complexity of the algorithm
is proportional to $N^{1.0691}2^{N/2}$ with $N$ being the size of the problem.
We also show that the recursive algorithm can be scaled up to HPC use cases
making feasible the simulation of threshold GBS up to $35-40$ photon clicks
without the needs of large-scale computational capacities.
