We present the non-perturbative computation of the entropy density in QCD for temperatures
ranging from 3 GeV up to the electro-weak scale, using $N_f=3$ flavours of massless O$(a)$-improved Wilson fermions.
We adopt a new strategy designed to be computationally efficient and based on formulating thermal QCD in a moving reference frame, where the
fields satisfy shifted boundary conditions in the temporal direction and periodic boundary conditions along the spatial ones.
In this setup the entropy density can be computed as the derivative of the free-energy density with respect to the shift parameter.
For each physical temperature, we perform Monte Carlo simulations at four values of the lattice spacing in order to extrapolate the
numerical data of the entropy density to the continuum limit. We achieve a final accuracy of approximatively $0.5$-$1.0\%$ and our results are
compared with predictions from high-temperature perturbation theory.