A fast C++ implementation of thermal functions

Solution method: We implement numerical quadrature, a Bessel function representation, asymptotic solutions expressed in terms of zeta functions and polylogarithms, and approximations and limits.