  Rapid and accurate evaluation of the nonlinear matter power spectrum, $P(k)$,
as a function of cosmological parameters and redshift is of fundamental
importance in cosmology. Analytic approximations provide an interpretable
solution, yet current approximations are neither fast nor accurate relative to
numerical emulators. We use symbolic regression to obtain simple analytic
approximations to the nonlinear scale, $k_\sigma$, the effective spectral
index, $n_{\rm eff}$, and the curvature, $C$, which are required for the
halofit model. We then re-optimise the coefficients of halofit to fit a wide
range of cosmologies and redshifts. We explore the space of analytic
expressions to fit the residuals between $P(k)$ and the optimised predictions
of halofit. Our results are designed to match the predictions of
EuclidEmulator2, but are validated against $N$-body simulations. Our symbolic
expressions for $k_\sigma$, $n_{\rm eff}$ and $C$ have root mean squared
fractional errors of 0.8%, 0.2% and 0.3%, respectively, for redshifts below 3
and a wide range of cosmologies. The re-optimised halofit parameters reduce the
root mean squared fractional error (compared to EuclidEmulator2) from 3% to
below 2% for wavenumbers $k=9\times10^{-3}-9 \, h{\rm Mpc^{-1}}$. We introduce
syren-halofit (symbolic-regression-enhanced halofit), an extension to halofit
containing a short symbolic correction which improves this error to 1%. Our
method is 2350 and 3170 times faster than current halofit and hmcode
implementations, respectively, and 2680 and 64 times faster than
EuclidEmulator2 (which requires running class) and the BACCO emulator. We
obtain comparable accuracy to EuclidEmulator2 and BACCO when tested on $N$-body
simulations. Our work greatly increases the speed and accuracy of symbolic
approximations to $P(k)$, making them significantly faster than their numerical
counterparts without loss of accuracy.
