We investigate the topological properties of N f = 2 + 1 QCD with physical quark masses, at temperatures around 500 MeV. With the aim of obtaining a reliable sampling of topological modes in a regime where the fluctuations of the topological charge Q are very rare, we adopt a multicanonical approach, adding a bias potential to the action which enhances the probability of suppressed topological sectors. This method permits to gain up to three orders magnitude in computational power in the explored temperature regime. Results at different lattice spacings and physical spatial volumes reveal no significant finite size effects and the presence, instead, of large finite cut-off effects, with the topological susceptibility which decreases by 3-4 orders of magnitude while moving from a ≃ 0.06 fm towards the continuum limit. The continuum extrapolation is in agreeement with previous lattice determinations with smaller uncertainties but obtained based on ansatzes justified by several theoretical assumptions. The parameter b 2 , related to the fourth order coefficient in the Taylor expansion of the free energy density f (θ), has instead a smooth continuum extrapolation which is in agreement with the dilute instanton gas approximation (DIGA); moreover, a direct measurement of the relative weights of the different topological sectors gives an even stronger support to the validity of DIGA.