  Perturbation theory (PT) is often used to model statistical observables
capturing the translation and rotation-invariant information in cosmological
density fields. PT produces higher-order corrections by integration over linear
statistics of the density fields weighted by kernels resulting from recursive
solution of the fluid equations. These integrals quickly become
high-dimensional and naively require increasing computational resources the
higher the order of the corrections. Here we show how to decouple the
integrands that often produce this issue, enabling PT corrections to be
computed as a sum of products of independent 1-D integrals. Our approach is
related to a commonly used method for calculating multi-loop Feynman integrals
in Quantum Field Theory, the Gegenbauer Polynomial $x$-Space Technique (GPxT).
We explicitly reduce the three terms entering the 2-loop power spectrum,
formally requiring 9-D integrations, to sums over successive 1-D radial
integrals. These 1-D integrals can further be performed as convolutions,
rendering the scaling of this method $N_{\rm g} \log N_{\rm g}$ with $N_{\rm
g}$ the number of grid points used for each Fast Fourier Transform. This method
should be highly enabling for upcoming large-scale structure redshift surveys
where model predictions at an enormous number of cosmological parameter
combinations will be required by Monte Carlo Markov Chain searches for the
best-fit values.
