Suppose $a>0$ and consider the integral $$ I(a) := \mathrm{p.v.}\int_0^\infty \frac{x}{(x^2 - a^2)( e^{x} - 1 )}\ dx $$ which is a principal value integral (there is a singularity at $x=a$).
How does one compute this? I thought about converting this to a different integral using $\int_0^\infty \cosh( a s ) e^{- s x} ds = \frac{x}{x^2- a^2}$, but this is only valid for $x>1$.
EDIT 1: Found an answer which I think is right (but my numerics don't confirm the answer, so some help would be appreciated): $$ I(a) = \frac{1}{2} \ln \left( \frac{a}{2\pi} \right) - \frac{1}{2} \mathrm{Re}\left( \psi \left( \frac{i a}{2\pi} \right) \right) $$ where $\psi$ is the polygamma-0 function. This does not match numerics from mathematica though, oddly enough. The way I derive this is to let $\epsilon > 0$ and use $$ \int_0^\infty \cos( a s) \sin( s x ) e^{- \epsilon s} ds = \frac{x (\epsilon^2 (a^2 + x^2 + \epsilon^2 ))}{( (x - a)^2 + \epsilon^2 )( (x+a)^2 + \epsilon^2 ) )} $$ which approaches $\frac{x}{x^2 - a^2}$ for $\epsilon \to 0$. $I(a)$ is now a double integral over $s$ and $x$, and the integrand falls off exponentially in every direction, so the modulus of the integrand converges. By Fubini's theorem, we can switch the order of integration. Using $\int_0^{\infty} \frac{\sin(s x)}{e^{x}-1} dx = \frac{\pi}{2} \mathrm{coth}\left( \pi s\right) - \frac{1}{2s}$ and then after that using $\int_0^\infty e^{- z y } [ s^{-1} - \mathrm{coth(s)}] = \psi( z/2 ) - \log(z/2)+1/z$ (taking the real part of that with $z= \frac{1}{\pi} ( a i + \epsilon )$, and then taking $\epsilon \to 0$) I get my answer.
EDIT 2: The $\psi$ function has an integral representation for any $z \in \mathbb{C}$ with $\mathrm{Re}[z]>0$ such that $$ \psi(z) = \log(z) - \frac{1}{2z} - 2 \int_0^\infty \frac{t \ dt}{(t^2 + z^2)(e^{2\pi t} - 1)} $$ which implies for any $w \in \mathbb{C}$ $$ \int_0^\infty \frac{x \ dx}{(x^2 + w^2)(e^{x} - 1)} = \frac{1}{2} \log\left(\frac{w}{2\pi}\right) - \frac{\pi}{2w} - \frac{1}{2} \psi\left(\frac{w}{2\pi}\right) \ . $$ Taking $w = \epsilon + i a $ for $\epsilon > 0$, and then taking the real part of each side in the limit $\epsilon \to 0^{+}$ seems to give the answer I quote above.