We study the deuteron production via the deuteron pion and nucleon catalysis reactions, πpn ↔ πd and N pn ↔ N d, by employing stochastic multi-particle reactions in the hadronic transport approach SMASH for the first time. This is an improvement compared to previous studies, which introduced an artificial fake resonance d to simulate these 3 ↔ 2 reactions as a chain of 2 ↔ 2 reactions. The derivation of the stochastic criterion for multi-particle reactions is presented in a comprehensive fashion and its implementation is tested against an analytic expression for the scattering rate and the equilibrating particle yields in box calculations. We then study Au + Au collisions at √ sNN = 7.7 GeV, where we find that multi-particle collisions substantially reduce the time required for deuterons to reach partial chemical equilibrium with nucleons. Subsequently, the final yield of d is practically independent from the number of d at particlization, confirming the results of previous studies. The mean transverse momentum and the integrated elliptic flow as a function of centrality are rather insensitive to the exact realization of the 2 ↔ 3 reactions.