We compute the three-loop beta functions of long-range multi-scalar models with general quartic interactions. The long-range nature of the models is encoded in a kinetic term with a Laplacian to the power $0<zeta<1$, rendering the computation of Feynman diagrams much harder than in the usual short-range case ($zeta=1$). As a consequence, previous results stopped at two loops, while six-loop results are available for short-range models. We push the renormalization group analysis to three loops, in an $epsilon=4zeta-d$ expansion at fixed dimension $d<4$, extensively using the Mellin-Barnes representation of Feynman amplitudes in the Schwinger parametrization. We then specialize the beta functions to various models with different symmetry groups: $O(N)$, $(mathbb{Z}_2)^N rtimes S_N$, and $O(N)times O(M)$. For such models, we compute the fixed points and critical exponents.