We present a precise computation of the topological charge distribution in the $SU(3)$ Yang-Mills theory. It is carried out on the lattice with high statistics Monte Carlo simulations by employing the clover discretization of the field strength tensor combined with the Yang-Mills gradient flow. The flow equations are integrated numerically by a fourth-order structure-preserving Runge-Kutta method. We have performed simulations at four lattice spacings and several lattice sizes to remove with confidence the systematic errors in the second (topological susceptibility $\chi_t^\text{YM}$) and the fourth cumulant of the distribution. In the continuum we obtain the preliminary results $t_0^2\chi_t^\text{YM}=6.53(8)\times 10^{-4}$ and the ratio between the fourth and the second cumulant $R=0.233(45)$. Our results disfavour the $\theta$-behaviour of the vacuum energy predicted by dilute instanton models, while they are compatible with the expectation from the large-$N_c$ expansion.