Integral $I=\int_0^\infty \frac{\ln(1+x)\ln(1+x^{-2})}{x} dx$

The infinite sum in Chen Wang's answer, that is, $ \displaystyle \sum_{n=1}^{\infty} \frac{H_{4n}}{n^{2}}$, can be evaluated using contour integration by considering the function $$f(z) = \frac{\pi \cot(\pi z) [\gamma + \psi(-4z)]}{z^{2}}, $$

where $\psi(z)$ is the digamma function and $\gamma$ is the Euler-Mascheroni constant.

The function $f(z)$ has poles of order $2$ at the positive integers, simple poles at the negative integers, simple poles at the positive quarter-integers, and a pole of order $4$ at the origin.

The function $\psi(-4z)$ does have simple poles at the positive half-integers, but they are cancelled by the zeros of $\cot( \pi z)$.

Now consider a square on the complex plane (call it $C_{N}$) with vertices at $\pm (N + \frac{1}{2}) \pm i (N +\frac{1}{2})$.

On the sides of the square, $\cot (\pi z)$ is uniformly bounded.

And when $z$ is large in magnitude and not on the positive real axis, $\psi(-4z) \sim \ln(-4z)$.

So $ \displaystyle \int_{C_{N}} f(z) \ dz $ vanishes as $N \to \infty$ through the positive integers.

Therefore,

$$\sum_{n=1}^{\infty} \text{Res} [f(z), n] + \sum_{n=1}^{\infty} \text{Res}[f(z),-n] + \text{Res}[f(z),0] + \sum_{n=0}^{\infty} \text{Res}\Big[f(z), \frac{2n+1}{4} \Big] =0 .$$

To determine the residues, we need the following Laurent expansions.

At the positive integers,

$$ \gamma + \psi (-4z) = \frac{1}{4} \frac{1}{z-n} + H_{4n} + \mathcal{O}(z-n) $$

and

$$ \pi \cot (\pi z) = \frac{1}{z-n} + \mathcal{O}(z-n) .$$

At the origin,

$$ \gamma+ \psi(-4z) = \frac{1}{4z} -4 \zeta(2) z -16 \zeta(3) z^{2} + \mathcal{O}(z^{3})$$

and $$ \pi \cot (\pi z) = \frac{1}{z} - 2 \zeta(2) z + \mathcal{O}(z^{3}) .$$

And at the positive quarter-integers,

$$ \gamma + \psi(-4z) = \frac{1}{4} \frac{1}{z-\frac{2n+1}{4}} + \mathcal{O}(1)$$

and

$$ \pi \cot (\pi z) = (-1)^{n} \pi + \mathcal{O}\Big(z- \frac{2n+1}{4} \Big) .$$

Then at the positive integers,

$$f(z) = \frac{1}{z^{2}} \Big( \frac{1}{4} \frac{1}{(z-n)^{2}} + \frac{H_{4n}}{z-n} + \mathcal{O}(1) \Big), $$

which implies

$$\begin{align} \text{Res} [f(z),n] &= \text{Res} \Big[ \frac{1}{4z^{2}} \frac{1}{(z-n)^{2}} , n \Big] + \text{Res} \Big[ \frac{1}{z^{2}} \frac{H_{4n}}{z-n}, n \Big] \\ &= - \frac{1}{2n^{3}} + \frac{H_{4n}}{n^{2}} .\end{align}$$

At the negative integers,

$$ \text{Res}[f(z),-n] = \frac{\gamma + \psi(4n)}{n^{2}} = \frac{H_{4n-1}}{n^{2}} = \frac{H_{4n}}{n^{2}} - \frac{1}{4n^{3}} . $$

At the origin,

$$ f(z) = \frac{1}{z^{2}} \Big( \frac{1}{4z^{2}} - \frac{\zeta(2)}{2} - 4 \zeta(2) - 16 \zeta(3) z + \mathcal{O}(z^{2}) \Big),$$

which implies

$$\text{Res}[f(z),0] = -16 \zeta(3) .$$

And at the positive quarter-integers,

$$ f(z) = \frac{\pi}{4z^{2}} \frac{(-1)^{n}}{z- \frac{2n+1}{4}} + \mathcal{O}(1),$$

which implies

$$ \begin{align} \text{Res} \Big[ f(z),\frac{2n+1}{4} \Big] &= \text{Res} \Big[\frac{\pi}{4z^{2}} \frac{(-1)^n}{z- \frac{2n+1}{4}}, \frac{2n+1}{4} \Big] \\ &= 4 \pi \ \frac{(-1)^{n}}{(2n+1)^{2}} . \end{align} $$

Putting everything together, we have

$$ - \frac{1}{2} \sum_{n=1}^{\infty} \frac{1}{n^{3}} + 2 \sum_{n=1}^{\infty} \frac{H_{4n}}{n^{2}} - \frac{1}{4} \sum_{n=1}^{\infty} \frac{1}{n^{3}} - 16 \zeta(3) + 4 \pi \sum_{n=0}^{\infty} \frac{(-1)^{n}}{(2n+1)^{2}} $$

$$ = - \frac{1}{2} \zeta(3) + 2 \sum_{n=1}^{\infty} \frac{H_{4n}}{n^{2}} - \frac{1}{4} \zeta(3) - 16 \zeta(3) + 4 \pi G = 0 .$$

Therefore,

$$ \sum_{n=1}^{\infty} \frac{H_{4n}}{n^{2}} = \frac{67}{8} \zeta(3) - 2 \pi G .$$

EDIT:

I found the Laurent expansion of $\psi(-4z)$ at the positive integers by using the functional equation of the digamma function to express $\psi(4z)$ as

$$ \psi(4z) = \psi(4z+4n+1) - \frac{1}{4z+4n} - \frac{1}{4z+4n-1} - \ldots - \frac{1}{4z} .$$

Then I evaluated the limit $$\lim_{z \to -n} (z+n) \psi(4z) = - \frac{1}{4}$$ and the limit $$\lim_{z \to -n} \Big(\psi(4z) + \frac{1}{4} \frac{1}{z+n} \Big) = - \gamma +H_{4n} .$$

This leads to the expansion $$\gamma + \psi (-4z) = \frac{1}{4} \frac{1}{z-n} + H_{4n} + \mathcal{O}(z-n) .$$

I did something similar to find the expansion at the positive quarter-integers.

I have a simple solution for this problem. Let $$ I(\alpha,\beta)=\int_0^\infty\frac{\ln(1+\alpha x)\ln(1+\beta x^{-2})}{x}dx. $$ Then $I(0,0)=I(\alpha,0)=I(0,\beta)=0$ and $I(1,1)=I$. It is easy to check \begin{eqnarray*} \frac{\partial^2 I}{\partial \alpha\partial\beta}&=&\int_0^\infty\frac{1}{(1+\alpha x)(x^2+\beta)}dx\\ &=&\frac{1}{2}\left(\frac{\pi}{\sqrt{\beta}}+2\alpha\ln\alpha+\alpha\ln\beta\right)\frac{1}{1+\alpha^2\beta}\\ &=&\frac{1}{2}\left(\frac{\pi}{\sqrt{\beta}}+2\alpha\ln\alpha+\alpha\ln\beta\right)\sum_{n=0}^\infty(-1)^n\alpha^{2n}\beta^n\\ &=&\frac{1}{2}\left(\pi\sum_{n=0}^\infty(-1)^n\alpha^{2n}\beta^{n-\frac{1}{2}}+2\sum_{n=0}^\infty(-1)^n(\alpha^{2n+1}\ln\alpha)\beta^n+\sum_{n=0}^\infty(-1)^n\alpha^{2n+1}\beta^n\ln\beta\right). \end{eqnarray*} Thus \begin{eqnarray*} I&=&\int_0^1\int_0^1\frac{\partial^2 I}{\partial \alpha\partial\beta}d\beta d\alpha\\ &=&\frac{1}{2}\int_0^1\left(\pi\sum_{n=0}^\infty\frac{(-1)^n}{n+\frac{1}{2}}\alpha^{2n}+2\sum_{n=0}^\infty\frac{(-1)^n}{n+1}\alpha^{2n+1}\ln\alpha-\sum_{n=0}^\infty\frac{(-1)^n}{(n+1)^2}\alpha^{2n+1}\right)d\alpha\\ &=&\pi\sum_{n=0}^\infty\frac{(-1)^n}{(2n+1)^2}-\sum_{n=0}^\infty\frac{(-1)^n}{4(n+1)^3}-\frac{1}{2}\sum_{n=0}^\infty\frac{(-1)^n}{2(n+1)^3}\\ &=&\pi G-\frac{1}{2}\sum_{n=0}^\infty\frac{(-1)^n}{(n+1)^3}\\ &=&\pi G-\frac{3\zeta(3)}{8}. \end{eqnarray*} Here we use $$\int_0^1 x^n\ln x dx=-\frac{1}{(n+1)^2}. $$

Incomplete answer:

$$ \begin{align*} I&=\int^{\infty}_{0}\frac{\log(1+x)\log(1+x^{-2})}{x}dx\\ &=\int^{1}_{0}\frac{\log(1+x)\log(1+x^{-2})}{x}dx+\int^{\infty}_{1}\frac{\log(1+x)\log(1+x^{-2})}{x}dx\\ &=\int^{1}_{0}\frac{\log(1+x)\log(1+x^{-2})}{x}dx+\int^{1}_{0}\frac{\log(1+x^{-1})\log(1+x^2)}{x}dx\\ &=2\underbrace{\int^{1}_{0}\frac{\log(1+x)\log(1+x^2)}{x}dx}_{=I_1}-\underbrace{\int^{1}_{0}\frac{\log(x)(2\log(1+x)+\log(1+x^2))}{x}dx}_{=I_2}.\\ \end{align*} $$

$$ \begin{align*} I_2&=\int^{1}_{0}\frac{\log(x)(2\log(1+x)+\log(1+x^2))}{x}dx\\ &=\sum^{\infty}_{n=1}\frac{(-1)^{n+1}}{n}\int^{1}_{0}\frac{\log(x)(2x^n+x^{2n})}{x}dx\\ &=\sum^{\infty}_{n=1}\frac{(-1)^{n+1}}{n}\left(-\frac{2}{n^2}-\frac{1}{(2n)^2}\right)\\ &=\frac{27}{16}\zeta(3). \end{align*} $$

$$ \begin{align*} I_1&=\int^{1}_{0}\frac{\log(1+x)\log(1+x^2)}{x}dx\\ &=\sum^{\infty}_{m,n=1}\frac{(-1)^{m+n}}{mn}\int^{1}_{0}x^{2m+n-1}dx\\ &=\sum^{\infty}_{m,n=1}\frac{(-1)^{m+n}}{mn(2m+n)}\\ &=\sum^{\infty}_{m=1}\frac{(-1)^{m}}{m}\sum^{\infty}_{n=1}\frac{(-1)^{n}}{n(2m+n)}\\ &=\sum^{\infty}_{m=1}\frac{(-1)^{m}}{2m^2}\sum^{\infty}_{n=1}(-1)^{n}\left(\frac{1}{n}-\frac{1}{2m+n}\right)\\ &=\sum^{\infty}_{m=1}\frac{(-1)^{m}}{2m^2}\left(\sum^{\infty}_{n=1}\frac{(-1)^{n}}{n}-\sum^{\infty}_{n=1}\frac{(-1)^{2m+n}}{2m+n}\right)\\ &=\sum^{\infty}_{m=1}\frac{(-1)^{m}}{2m^2}\sum^{2m}_{n=1}\frac{(-1)^{n}}{n}\\ &=\sum^{\infty}_{m=1}\frac{(-1)^{m-1}(H_{2m}-H_{m})}{2m^2}\\ &=\sum^{\infty}_{m=1}\frac{H_{2m}-H_{m}}{2m^2}-2\sum^{\infty}_{m=1}\frac{H_{4m}-H_{2m}}{2(2m)^2}\\ &=\frac14\sum^{\infty}_{m=1}\frac{-H_{4m}+3H_{2m}-2H_{m}}{m^2}\\ &=-\frac14\sum^{\infty}_{m=1}\frac{H_{4m}}{m^2}+\frac{17}{16}\zeta(3) \end{align*}$$

Therefore, $I=2I_1-I_2=-\frac12\sum^{\infty}_{m=1}\frac{H_{4m}}{m^2}+\frac{61}{16}\zeta(3)$.

It remains to prove that $\sum^{\infty}_{m=1}\frac{H_{4m}}{m^2}\stackrel?=\frac{67}{8}\zeta(3)-2\pi G$, but I haven't found a way yet.