I seem to have found a more "interesting," way to compute the integral in the earlier post (not the general form though).
$$I = \int_{0}^{1} \log(x)\log(1-x) dx$$
We can use:
$$\log(1-x) = -\sum_{n=1}^{\infty} \frac{x^n}{n}$$
$$\log(x)\log(1-x) = -\sum_{n=1}^{\infty} \frac{\log(x)x^n}{n}$$
$$\int_{0}^{1} \log(x)\log(1-x) = -\sum_{n=1}^{\infty} \frac{1}{n}\cdot \int_{0}^{1} \log(x) x^n dx$$
With
differentiation under the integral sign (leibniz rule) trick for the integral in the RHS (with differentiation $d/dn x^n = \log(x)x^n$), it is easy to see that:
$$\int_{0}^{1} \log(x)\log(1-x) dx = \sum_{n=1}^{\infty} \frac{1}{n(n+1)^2} = I$$
First consider a square contour $C$ below:
View attachment 3848
I am assuming, as $R \to \infty$ we have: $\displaystyle \oint_{C} f(z) dz \to 0$.
I do not have a proof of this, any proof is welcome, but I will work with this assumption.Consider:
$$f(z) = \frac{\psi(-z)}{z(1+z)^2}$$
$$\mathrm{Res}_{z=n} f(z) = \mathrm{Res} \psi(-z) \cdot \frac{1}{n(n+1)^2} = \frac{1}{n(n+1)^2}$$
You can check that the residue of $\psi(-z) = 1$ by using series and finding the coefficient of $\frac{1}{z-n}$
Because $z=n$ for $n \ge 1$ we only have a simple pole (order 1).
$$\mathrm{Res}_{z=0} f(z) = -2 - \gamma$$ we show this by:
Digamma function - Wikipedia, the free encyclopedia
$$\psi(-z) = \frac{1}{z} - \gamma - \zeta(2)z + ...$$
$$\frac{1}{z(1+z)^2} = \frac{1}{z} - 2 + 3x + ...$$
Multiplying the two series we receive:
$$\frac{\psi(-z)}{z(1+z)^2} = \left( \frac{1}{z} - \gamma - \zeta(2)z + ... \right) \cdot \left( \frac{1}{z} - 2 + 3z + ... \right)$$
$$ = \left( \frac{1}{z^2} -\frac{2}{z} + 3 + ... \right) + \left( \frac{1}{z^2} -\frac{\gamma}{z}+ 2\zeta(2) + ... \right)$$
$$ = ... + \frac{1}{z}\cdot \left( -\gamma - 2 \right) + ...$$
Hence,
$$\mathrm{Res}_{z=0} f(z) = -\gamma - 2$$
Next we find:
$$\mathrm{Res}_{z=-1} f(z) = \gamma + \zeta(2)$$
This is a double pole we simply use the formula:
$$\mathrm{Res}_{z=-1} f(z) = \lim_{z = -1} - \frac{\psi(-z) + z\psi_{1}(-z)}{z^2} = -\psi(1) + \psi_{1}(1) = \gamma + \zeta(2)$$
The values for $\psi(1)$ and $\psi_{1}(1)$ I took from Wikipedia Tables, but I suppose ti s derivable by the series representations...
Anyway.
$$\frac{1}{2\pi i} \cdot \oint_{C} f(z) dz = \sum \text{Residues in contour} = 0$$
$$\sum_{n=1}^{\infty} \frac{1}{n(n+1)^2} -(\gamma + 2) + \gamma + \zeta(2) = 0$$
$$\sum_{n=1}^{\infty} \frac{1}{n(n+1)^2} = 2 - \zeta(2)$$
Finally, just to make sure of a clean ending:
$$I = \int_{0}^{1} \log(x)\log(1-x) dx = \sum_{n=1}^{\infty} \frac{1}{n(n+1)^2} = 2 - \zeta(2) \space \space \space \space \blacksquare$$