We introduce a method to numerically compute equilibrium measures for
problems with attractive-repulsive power law kernels of the form $K(x-y) =
\frac{|x-y|^\alpha}{\alpha}-\frac{|x-y|^\beta}{\beta}$ using recursively
generated banded and approximately banded operators acting on expansions in
ultraspherical polynomial bases. The proposed method reduces what is naively a
difficult to approach optimization problem over a measure space to a
straightforward optimization problem over one or two variables fixing the
support of the equilibrium measure. The structure and rapid convergence
properties of the obtained operators results in high computational efficiency
in the individual optimization steps. We discuss stability and convergence of
the method under a Tikhonov regularization and use an implementation to
showcase comparisons with analytically known solutions as well as discrete
particle simulations. Finally, we numerically explore open questions with
respect to existence and uniqueness of equilibrium measures as well as gap
forming behaviour in parameter ranges of interest for power law kernels, where
the support of the equilibrium measure splits into two intervals.