  Calculating fusion reactivity involves a complex six-dimensional integral of
the fusion cross section and ion velocity distributions of two reactants. We
demonstrate a simple Monte Carlo approach that efficiently computes this
integral for arbitrary ion velocity distributions with a time complexity of
$O(N)$, where $N$ is the number of samples. This approach generates random
numbers that satisfy the reactant velocity distributions. In cases where these
numbers are not readily available, we propose using Gaussian random numbers
with weighted factors. For cases where only a small number of $N$ samples are
available, a $O(N^2)$ method can be used. We benchmarked this approach against
analytical results for drift bi-Maxwellian distributions and provided examples
of drift ring beam and slowing down distributions. Our results show that the
error can be less than 1\% with $N\sim10^4$ samples for our standard approach.
