I have been trying to evaluate the following family of integrals:
f:(0,\infty)^2 \rightarrow \mathbb{R} \, , \, f(\alpha,\beta) = \int \limits_0^\infty \frac{\ln (1+x^\alpha) \ln (1+x^{-\beta})}{x} \, \mathrm{d} x \, .
The changes of variables \frac{1}{x} \rightarrow x, x^\alpha \rightarrow x and x^\beta \rightarrow x yield the symmetry properties
\tag{1} f(\alpha,\beta) = f(\beta,\alpha) = \frac{1}{\alpha} f\left(1,\frac{\beta}{\alpha}\right) = \frac{1}{\alpha} f\left(\frac{\beta}{\alpha},1\right) = \frac{1}{\beta} f\left(\frac{\alpha}{\beta},1\right) = \frac{1}{\beta} f\left(1,\frac{\alpha}{\beta}\right)
for \alpha,\beta > 0 .
Using this result one readily computes f(1,1) = 2 \zeta (3) . Then (1) implies that
f(\alpha,\alpha) = \frac{2}{\alpha} \zeta (3)
holds for \alpha > 0 . Every other case can be reduced to finding f(1,\gamma) for \gamma > 1 using (1).
An approach based on xpaul's answer to this question employs Tonelli's theorem to write
\tag{2} f(1, \gamma) = \int \limits_0^\infty \int \limits_0^1 \int \limits_0^1 \frac{\mathrm{d}u \, \mathrm{d}v \, \mathrm{d}x}{(1+ux)(v+x^\gamma)} = \int \limits_0^1 \int \limits_0^1 \int \limits_0^\infty \frac{\mathrm{d}x \, \mathrm{d}u \, \mathrm{d}v}{(1+ux)(v+x^\gamma)} \, .
The special case f(1,2) = \pi \mathrm{C} - \frac{3}{8} \zeta (3) is then derived via partial fraction decomposition (\mathrm{C} is Catalan's constant). This technique should work at least for \gamma \in \mathbb{N} (it also provides an alternative way to find f(1,1)), but I would imagine that the calculations become increasingly complicated for larger \gamma .
Mathematica manages to evaluate f(1,\gamma) in terms of \mathrm{C}, \zeta(3) and an acceptably nice finite sum of values of the trigamma function \psi_1 for some small, rational values of \gamma > 1 (before resorting to expressions involving the Meijer G-function for larger arguments). This gives me some hope for a general formula, though I have not yet been able to recognise a pattern.
Therefore my question is:
How can we compute f(1,\gamma) for general (or at least integer/rational) values of \gamma > 1 ?
Update 1:
Symbolic and numerical evaluations with Mathematica strongly suggest that
f(1, n) = \frac{1}{n (2 \pi)^{n-1}} \mathrm{G}_{n+3, n+3}^{n+3,n+1} \left(\begin{matrix} 0, 0, \frac{1}{n}, \dots, \frac{n-1}{n}, 1 , 1 \\ 0,0,0,0,\frac{1}{n}, \dots, \frac{n-1}{n} \end{matrix} \middle| \, 1 \right)
holds for n \in \mathbb{N} . These values of the Meijer G-function admit an evaluation in terms of \zeta(3) and \psi_1 \left(\frac{1}{n}\right), \dots, \psi_1 \left(\frac{n-1}{n}\right) at least for small (but likely all) n \in \mathbb{N} .
Interesting side note: The limit
\lim_{\gamma \rightarrow \infty} f(1,\gamma+1) - f(1,\gamma) = \frac{3}{4} \zeta(3)
follows from the definition.
Update 2:
Assume that m, n \in \mathbb{N} are relatively prime (i.e. \gcd(m,n) = 1). Then the expression for f(m,n) given in Sangchul Lee's answer can be reduced to
\begin{align} f(m,n) &= \frac{2}{m^2 n^2} \operatorname{Li}_3 ((-1)^{m+n}) \\ &\phantom{=} - \frac{\pi}{4 m^2 n} \sum \limits_{j=1}^{m-1} (-1)^j \csc\left(j \frac{n}{m} \pi \right) \left[\psi_1 \left(\frac{j}{2m}\right) + (-1)^{m+n} \psi_1 \left(\frac{m + j}{2m}\right) \right] \\ &\phantom{=} - \frac{\pi}{4 n^2 m} \sum \limits_{k=1}^{n-1} (-1)^k \csc\left(k \frac{m}{n} \pi \right) \left[\psi_1 \left(\frac{k}{2n}\right) + (-1)^{n+m} \psi_1 \left(\frac{n + k}{2n}\right) \right] \\ &\equiv F(m,n) \, . \end{align}
Further simplifications depend on the parity of m and n.
This result can be used to obtain a solution for arbitrary rational arguments: For \frac{n_1}{d_1} , \frac{n_2}{d_2} \in \mathbb{Q}^+ equation (1) yields
\begin{align} f\left(\frac{n_1}{d_1},\frac{n_2}{d_2}\right) &= \frac{d_1}{n_1} f \left(1,\frac{n_2 d_1}{n_1 d_2}\right) = \frac{d_1}{n_1} f \left(1,\frac{n_2 d_1 / \gcd(n_1 d_2,n_2 d_1)}{n_1 d_2 / \gcd(n_1 d_2,n_2 d_1)}\right) \\ &= \frac{d_1 d_2}{\gcd(n_1 d_2,n_2 d_1)} f\left(\frac{n_1 d_2}{\gcd(n_1 d_2,n_2 d_1)},\frac{n_2 d_1}{\gcd(n_1 d_2,n_2 d_1)}\right) \\ &= \frac{d_1 d_2}{\gcd(n_1 d_2,n_2 d_1)} F\left(\frac{n_1 d_2}{\gcd(n_1 d_2,n_2 d_1)},\frac{n_2 d_1}{\gcd(n_1 d_2,n_2 d_1)}\right) \, . \end{align}
Therefore I consider the problem solved in the case of rational arguments. Irrational arguments can be approximated by fractions, but if anyone can come up with a general solution: you are most welcome to share it. ;)
Only a comment. We have
\int_{0}^{\infty} \frac{\log(1+\alpha x)\log(1+\beta/x)}{x} \, dx = 2\operatorname{Li}_3(\alpha\beta) - \operatorname{Li}_2(\alpha\beta)\log(\alpha\beta)
which is valid initially for \alpha, \beta > 0 and extends to a larger domain by the principle of analytic continuation. Then for integers m, n \geq 1 we obtain
\begin{align*} f(m, n) &=\int_{0}^{\infty} \frac{\log(1+x^m)\log(1+x^{-n})}{x}\,dx \\ &\hspace{6em} = \sum_{j=0}^{m-1}\sum_{k=0}^{n-1} \left[ 2\operatorname{Li}_3\left(e^{i(\alpha_j+\beta_k)}\right) - i(\alpha_j+\beta_k)\operatorname{Li}_2\left(e^{i(\alpha_j+\beta_k)}\right) \right], \end{align*}
where \alpha_j = \frac{2j-m+1}{n}\pi and \beta_k = \frac{2k-n+1}{n}\pi. (Although we cannot always split complex logarithms, this happens to work in the above situation.) By the multiplication formula, this simplifies to
\begin{align*} f(m, n) &= \frac{2\gcd(m,n)^3}{m^2n^2}\operatorname{Li}_3\left((-1)^{(m+n)/\gcd(m,n)}\right) \\ &\hspace{2em} - \frac{i}{n} \sum_{j=0}^{m-1} \alpha_j \operatorname{Li}_2\left((-1)^{n-1}e^{in\alpha_j}\right) \\ &\hspace{2em} - \frac{i}{m} \sum_{k=0}^{n-1} \beta_k \operatorname{Li}_2\left((-1)^{m-1}e^{im\beta_k}\right). \end{align*}
Here, \gcd(m,n) is the greatest common divisor of m and n.
The following code tests the above formula.
