A closed-form expression for $\int_0^\infty \frac{\ln (1+x^\alpha) \ln (1+x^{-\beta})}{x} \, \mathrm{d} x$ 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. ;)
 A: Only a comment.
For $\,z> 0\,$ I’ve got:

$$f(1,z) = \frac{3}{4}\zeta(3)\left(z+\frac{1}{z^2}\right) + 2 g(z)$$

with $\enspace\displaystyle g(z):=\sum\limits_{n=1}^\infty \frac{(-1)^{n-1}}{n} \sum\limits_{k=1}^\infty \frac{(-1)^{k-1}}{k(zn+k)}\enspace$ and $\enspace\displaystyle  g(\frac{1}{z})=zg(z)$  
Maybe someone can take it from here to create formulas.
Note: $\enspace\displaystyle  f(1,\frac{1}{z})=zf(1,z)\enspace$ which is equivalent to $\enspace\displaystyle \frac{1}{\alpha}f(1,\frac{\beta}{\alpha})=\frac{1}{\beta}f(1,\frac{\alpha}{\beta})$

Hint:
$\displaystyle zg(z)=\sum\limits_{n=1}^\infty \frac{(-1)^{n-1}}{n^2}\left(H(zn)- H(\frac{zn}{2})\right)=\int\limits_0^1\frac{Li_2(-t^z) - Li_2(-t^{z/2})}{1-t}dt $ 
$\hspace{1.1cm}\displaystyle = \frac{\pi^2\ln 2}{12} + \int\limits_0^1\frac{Li_2(-x^z)}{1+x}dx = z \int\limits_0^1\frac{\ln(1+x)\ln(1+x^z)}{x}$ 
with $\enspace\displaystyle H(x):=\sum\limits_{k=1}^\infty \frac{x}{k(x+k)}=\gamma + \psi(1+x) =\int\limits_0^1\frac{1-t^x}{1-t}dt$
It's trivial that $\,f(1,z)\,$ can be split:
$\hspace{1.1cm}\displaystyle \int\limits_0^1 \frac{\ln(1+x)\ln(1+x^{-z})}{x}dx = \frac{3z}{4}\zeta(3) + g(z) = g(-z)$
$\hspace{1.1cm}\displaystyle \int\limits_1^\infty \frac{\ln(1+x)\ln(1+x^{-z})}{x}dx = \frac{3}{4z^2}\zeta(3) + g(z)$
A: 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.

A: The integral of $x^p$ times a product of two Meijer G-functions gives a Fox H-function:
$$\int_0^\infty \frac {\ln(1+x^\alpha) \ln(1+x^{-\beta})} x dx =
\frac 1 \alpha \int_0^\infty
 G_{2,2}^{1,2} \left( x \middle| {1, 1 \atop 1, 0} \right)
 G_{2,2}^{2,1} \left( x^{\beta/\alpha} \middle| {0, 1 \atop 0, 0} \right) 
 \frac {dx} x = \\
\frac 1 \alpha H_{4,4}^{4,2} \left( 1 \middle|
 {(0, 1), (0, \frac \beta \alpha), (1, \frac \beta \alpha), (1, 1) \atop
  (0, 1), (0, 1), (0, \frac \beta \alpha), (0, \frac \beta \alpha)} \right) =
\frac 1 {2 \pi i } \int_\mathcal L
 \frac {\pi^2 \csc \pi s \csc \frac {\pi \beta s} \alpha} {\beta s^2} ds.$$
When $\beta/\alpha$ is rational, the H-function reduces to a G-function.
A: This is not an answer since just based on a CAS.
For integer values of $\gamma$, we have
$$f(1,1)=2 \zeta (3)$$
$$f(1,2)=\pi  C-\frac{3 \zeta (3)}{8}$$
$$f(1,3)=\frac{2}{27} \left(3 \zeta (3)+\sqrt{3} \pi  \left(\psi
   ^{(1)}\left(\frac{1}{3}\right)-\psi ^{(1)}\left(\frac{2}{3}\right)\right)\right)$$
$$f(1,4)=\frac{1}{64} \left(-16 \pi  C-6 \zeta (3)+\sqrt{2} \pi  \left(\psi
   ^{(1)}\left(\frac{1}{8}\right)+\psi ^{(1)}\left(\frac{3}{8}\right)-\psi
   ^{(1)}\left(\frac{5}{8}\right)-\psi ^{(1)}\left(\frac{7}{8}\right)\right)\right)$$
I have been able to produce the results up to $f(1,6)$ without noticing any pattern; for $\gamma > 6$ start to appear  Meijer G functions (which are not the most pleasant - at least to me).
For some rational values of $\gamma$, I (say the CAS) obtained
$$f\left(1,\frac{1}{2}\right)=2 \pi  C-\frac{3 \zeta (3)}{4}$$
For $f\left(1,\frac{3}{2}\right)$ (no room to put it on a line) 
$$\frac{1}{108} \left(-72 \pi  C-9 \zeta (3)+2 \sqrt{3} \pi  \left(\psi
   ^{(1)}\left(\frac{1}{6}\right)+\psi ^{(1)}\left(\frac{1}{3}\right)-\psi
   ^{(1)}\left(\frac{2}{3}\right)-\psi ^{(1)}\left(\frac{5}{6}\right)\right)\right)$$
For $f\left(1,\frac{5}{2}\right)$, I also got a result (quite long !).
$$f\left(1,\frac{1}{3}\right)=\frac{2}{9} \left(3 \zeta (3)+\sqrt{3} \pi  \left(\psi
   ^{(1)}\left(\frac{1}{3}\right)-\psi ^{(1)}\left(\frac{2}{3}\right)\right)\right)$$
Now, I really need a break (my computer too).
I hope and wish that this could be of some interest for you.
A: I have found a way to evaluate $f(1,n)$ for $n \in \mathbb{N}$ . Here's a sketch:
Start by letting $x = v^{1/n} t$ in equation $(2)$. Then
$$ f(1,n) = \int \limits_0^1 \int \limits_0^1 \int \limits_0^\infty \frac{\mathrm{d} t}{(t + (u v^{1/n})^{-1}) (1 + t^n)} \, \frac{\mathrm{d} u}{u} \, \frac{\mathrm{d} v}{v} \equiv \int \limits_0^1 \int \limits_0^1 X((u v^{1/n})^{-1}) \, \frac{\mathrm{d} u}{u} \, \frac{\mathrm{d} v}{v} \, .$$
A partial fraction decomposition yields
\begin{align}
X(w) &\equiv \int \limits_0^\infty \frac{\mathrm{d} t}{(t + w) (1 + t^n)} = \frac{1}{1+(-w)^n} \int \limits_0^\infty \left[ \frac{1}{t+w} - \frac{1}{1+t^n} \sum \limits_{k=0}^{n-1} (-w)^{n-1-k} t^k\right] \, \mathrm{d} t \\
&\equiv Y(w) + Z(w) \, .
\end{align}
Here
$$ Y(w) = \frac{1}{1+(-w)^n} \int \limits_0^\infty \left[ \frac{1}{t+w} - \frac{t^{n-1}}{1+t^n} \right] \, \mathrm{d} t = \frac{-\ln(w)}{1+(-w)^n} $$
and
$$ Z(w) = - \frac{1}{1+(-w)^n} \sum \limits_{k=0}^{n-2} (-w)^{n-1-k} \int \limits_0^\infty \frac{t^k}{1 + t^n} \, \mathrm{d} t = - \frac{\pi}{n} \sum \limits_{m=1}^{n-1} \csc\left(\frac{m}{n} \pi\right)\frac{(-w)^{n-m}}{1+(-w)^n} $$
have been introduced (the last step follows from this question). Now check that
\begin{align} 
\int \limits_0^1 \int \limits_0^1 Y((u v^{1/n})^{-1}) \, \frac{\mathrm{d} u}{u} \, \frac{\mathrm{d} v}{v} &= \frac{1}{n} \int \limits_0^1 \int \limits_0^1 \frac{\ln(u^n v) u^n v}{u^n v +(-1)^n} \, \frac{\mathrm{d} v}{v} \, \frac{\mathrm{d} u}{u} \\
&\stackrel{(\mathrm{a})}{=} \frac{(-1)^{n-1}}{n} \int \limits_0^1 \int \limits_{-n \ln(u)}^\infty \frac{s}{\mathrm{e}^s +(-1)^n} \, \mathrm{d} s \, \frac{\mathrm{d} u}{u} \\
&\stackrel{(\mathrm{b})}{=} \frac{(-1)^{n-1}}{n} \int \limits_0^\infty \frac{s}{\mathrm{e}^s +(-1)^n} \int \limits_{\mathrm{e}^{-s/n}}^1 \frac{\mathrm{d} u}{u} \, \mathrm{d} s \\
&= \frac{(-1)^{n-1}}{n^2} \int \limits_0^\infty \frac{s^2}{\mathrm{e}^s +(-1)^n} \, \mathrm{d} s \\
&\stackrel{(\mathrm{c})}{=} \begin{cases} \frac{2}{n^2} \zeta (3) \, , & n ~  \mathrm{odd} \\ - \frac{3}{2 n^2} \zeta (3) \, , & n ~  \mathrm{even} \end{cases}
\end{align}
holds (explanation: $(\mathrm{a})$ $v = u^{-n} \mathrm{e}^{-s}$ , $(\mathrm{b})$ Tonelli, $(\mathrm{c})$ integral representations of $\zeta$). Similarly,
\begin{align}
\int \limits_0^1 \int \limits_0^1 Z((u v^{1/n})^{-1}) \, \frac{\mathrm{d} u}{u} \, \frac{\mathrm{d} v}{v} &= - \frac{\pi}{n} \sum \limits_{m=1}^{n-1} (-1)^{n-m} \csc\left(\frac{m}{n} \pi\right) \int \limits_0^1 \int \limits_0^1 \frac{(u^n v)^{m/n}}{(-1)^n + u^n v} \, \frac{\mathrm{d} v}{v} \, \frac{\mathrm{d} u}{u} \\
&\stackrel{(\mathrm{d})}{=} - \frac{\pi}{n} \sum \limits_{m=1}^{n-1} (-1)^{n-m} \csc\left(\frac{m}{n} \pi\right) \int \limits_0^1 \int \limits_0^{u^n} \frac{s^{\frac{m}{n} -1}}{(-1)^n + s} \, \mathrm{d} s \, \frac{\mathrm{d} u}{u} \\
&\stackrel{(\mathrm{e})}{=} - \frac{\pi}{n^2} \sum \limits_{m=1}^{n-1} (-1)^{n-m} \csc\left(\frac{m}{n} \pi\right) \int \limits_0^1 \frac{- \ln(s) s^{\frac{m}{n} -1}}{(-1)^n + s} \, \mathrm{d} s  \\
&\stackrel{(\mathrm{f})}{=} - \frac{\pi}{n^2} \sum \limits_{m=1}^{n-1} (-1)^{n-m} \csc\left(\frac{m}{n} \pi\right) \begin{cases} - \operatorname{\psi_1} \left(\frac{m}{n}\right) , & \!\!\! n ~ \mathrm{odd} \\ \frac{1}{4} \left[\operatorname{\psi_1} \left(\frac{m}{2 n}\right) - \operatorname{\psi_1} \left(\frac{n+m}{2 n}\right) \right]  , & \!\!\! n ~ \mathrm{even} \end{cases} \\
\end{align}
can be computed (explanation: $(\mathrm{d})$ $v = u^{-n} s$ , $(\mathrm{e})$ Tonelli, $(\mathrm{f})$ integral representation for $\psi_1$).
Combine these results to find
$$ f(1,n) = \frac{2}{n^2} \zeta (3) - \frac{\pi}{n^2} \sum \limits_{m=1}^{n-1} (-1)^{m} \csc\left(\frac{m}{n} \pi\right) \operatorname{\psi_1} \left(\frac{m}{n}\right) $$
for odd and 
$$ f(1,n) = - \frac{3}{2 n^2} \zeta (3) - \frac{\pi}{4 n^2} \sum \limits_{m=1}^{n-1} (-1)^{m} \csc\left(\frac{m}{n} \pi\right) \left[\operatorname{\psi_1} \left(\frac{m}{2 n}\right) - \operatorname{\psi_1} \left(\frac{n+m}{2 n}\right) \right] $$
for even $n \in \mathbb{N}$ .
Judging from the result for $f(1,\frac{3}{2})$ in Claude Leibovici's answer, a similar approach should work for rational values of $\gamma$, but I am not that optimistic about irrational arguments.
