  Bosonic exchange symmetry leads to fascinating quantum phenomena, from
exciton condensation in quantum materials to the superfluidity of liquid
Helium-4. Unfortunately, path integral molecular dynamics (PIMD) simulations of
bosons are computationally prohibitive beyond $\mathord{\sim} 100$ particles,
due to a cubic scaling with the system size. We present an algorithm that
reduces the complexity from cubic to quadratic, allowing the first simulations
of thousands of bosons using PIMD. Our method is orders of magnitude faster,
with a speedup that scales linearly with the number of particles and the number
of imaginary time slices (beads). Simulations that would have otherwise taken
decades can now be done in days. In practice, the new algorithm eliminates most
of the added computational cost of including bosonic exchange effects, making
them almost as accessible as PIMD simulations of distinguishable particles.
