We report on the perturbative computation of the 3d lattice Yang-Mills free energy to four loops by means of Numerical Stochastic Perturbation Theory. The known first and second orders have been correctly reproduced; the third and fourth order coefficients are new results and the known logarithmic IR divergence in the fourth order has been correctly identified. Progress is being made in switching to the gluon mass IR regularization and the related inclusion of the Faddeev-Popov determinant.