  Computing the matter power spectrum, $P(k)$, as a function of cosmological
parameters can be prohibitively slow in cosmological analyses, hence emulating
this calculation is desirable. Previous analytic approximations are
insufficiently accurate for modern applications, so black-box, uninterpretable
emulators are often used. We utilise an efficient genetic programming based
symbolic regression framework to explore the space of potential mathematical
expressions which can approximate the power spectrum and $\sigma_8$. We learn
the ratio between an existing low-accuracy fitting function for $P(k)$ and that
obtained by solving the Boltzmann equations and thus still incorporate the
physics which motivated this earlier approximation. We obtain an analytic
approximation to the linear power spectrum with a root mean squared fractional
error of 0.2% between $k = 9\times10^{-3} - 9 \, h{\rm \, Mpc^{-1}}$ and across
a wide range of cosmological parameters, and we provide physical
interpretations for various terms in the expression. Our analytic approximation
is 950 times faster to evaluate than camb and 36 times faster than the neural
network based matter power spectrum emulator BACCO. We also provide a simple
analytic approximation for $\sigma_8$ with a similar accuracy, with a root mean
squared fractional error of just 0.1% when evaluated across the same range of
cosmologies. This function is easily invertible to obtain $A_{\rm s}$ as a
function of $\sigma_8$ and the other cosmological parameters, if preferred. It
is possible to obtain symbolic approximations to a seemingly complex function
at a precision required for current and future cosmological analyses without
resorting to deep-learning techniques, thus avoiding their black-box nature and
large number of parameters. Our emulator will be usable long after the codes on
which numerical approximations are built become outdated.
