We approximate boson thermodynamic integrals as polynomials in two variables chosen to give the correct limiting expansions and to smoothly interpolate into other regimes. With 10 free parameters, an accuracy of better than 0.009\% is achieved for the pressure, internal energy density and the number density. We also revisit the fermion case, originally addressed by Eggleton, Faulkner and Flannery (1973), and substantially improve the accuracy of their fits.