We derive the finite-temperature equation of state of dark matter superfluids with 2-body and 3-body contact interactions. The latter case is relevant to a recently proposed model of dark matter superfluidity that unifies the collisionless aspects of dark matter with the empirical success of MOdified Newtonian Dynamics at fitting galactic rotation curves. The calculation uses a self-consistent mean-field approximation. It relies on the Hartree-Fock-Bogoliubov approximation and follows the Yukalov-Yukalova proposal to circumvent the well-known Hohenberg-Martin dilemma. The resulting equation of state is consistent with a gapless spectrum, and simultaneously satisfies the equation of motion for the condensate wavefunction. As an application, we derive the finite-temperature density profile for dark matter superfluids, assuming spherical symmetry and uniform temperature. The density profiles consist of a nearly homogeneous superfluid core, surrounded by an isothermal "atmosphere" of normal particles, with the transition taking place at the critical density.