The Debye functions are defined like so: ${D_n\left(x\right)} = {\frac{n}{x^n} \cdot {\int_0^x{\frac{t^n}{e^t - 1}dt}}}$.
I'm trying to evaluate the functions for $n$ from one to four and for $\left\lvert x \right\rvert < 10$. The required number of evaluations is great, and I require hundreds of bits of accuracy.
Using Julia and the Julia package QuadGK, implementing the integral in the definition is simple enough, however it's slow for tiny values of $x$, and the integrand even overflows (due to the MPFR exponent not being wide enough) for very tiny values of $x$.
I wonder if there's a way compute $D_1$, $D_2$, $D_3$, $D_4$ without resorting to generic numeric integration, perhaps it would be possible to design a custom algorithm for evaluating the functions that would be faster?
expm1
instead of computingexp(t)-1
. $\endgroup$expm1
seems to solve all my problems so far. $\endgroup$