The surface free energies of seven different facets of tungsten (W) are obtained up to the melting point with full account of all the relevant thermal excitations; in particular, thermal atomic vibrations, electronic excitations, and their mutual coupling. The latter is done using ab initio molecular dynamics simulations coupled with the thermodynamic integration technique. In this way, the calculations contain almost no error but the one related to the used exchange-correlation functional, which makes the results truly first principles. The obtained results are compared with previous quasiharmonic calculations for the surface free energies of W and experimental data. The anharmonic contribution is, as expected, important for open surfaces at high temperatures, which leads to a temperature dependence of the surface energy anisotropy. The calculated Wulff shapes and surface energies are in excellent agreement with experimental data close to the melting point, where the crystalline structure of the surface layers is destroyed by a dramatic mobility of the atoms there.