  We present efficient algorithms for computing the $N$-point correlation
functions (NPCFs) of random fields in arbitrary $D$-dimensional homogeneous and
isotropic spaces. Such statistics appear throughout the physical sciences, and
provide a natural tool to describe stochastic processes. algorithms for
computing the NPCF components have $\mathcal{O}(n^N)$ complexity (for a data
set containing $n$ particles); their application is thus computationally
infeasible unless $N$ is small. By projecting the statistic onto a
suitably-defined angular basis, we show that the estimators can be written in a
separable form, with complexity $\mathcal{O}(n^2)$, or $\mathcal{O}(n_{\rm
g}\log n_{\rm g})$ if evaluated using a Fast Fourier Transform on a grid of
size $n_{\rm g}$. Our decomposition is built upon the $D$-dimensional
hyperspherical harmonics; these form a complete basis on the $(D-1)$-sphere and
are intrinsically related to angular momentum operators. Concatenation of
$(N-1)$ such harmonics gives states of definite combined angular momentum,
forming a natural separable basis for the NPCF. As $N$ and $D$ grow, the number
of basis components quickly becomes large, providing a practical limitation to
this (and all other) approaches: however, the dimensionality is greatly reduced
in the presence of symmetries; for example, isotropic correlation functions
require only states of zero combined angular momentum. We provide a
\textsc{Julia} package implementing our estimators, and show how they can be
applied to a variety of scenarios within cosmology and fluid dynamics. The
efficiency of such estimators will allow higher-order correlators to become a
standard tool in the analysis of random fields.
