Numerical approximation to the thermodynamic integrals

We approximate boson thermodynamic integrals as polynomials in two variables chosen to give the correct limiting expansion 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 number density. We also revisit the fermion case, originally addressed by Eggleton, Faulkner, & Flannery (1973), and substantially improve the accuracy of their fits.