Evaluating the integral $\int_1^{\infty} \int_1^{\infty} \frac{\Gamma(x+1)\Gamma(y+1)}{x y \Gamma(x+y+2)} \ dx \ dy$ During the study of some integrals I came across a very interesting integral, that is 
$$\int_1^{\infty} \int_1^{\infty} \frac{\Gamma(x+1)\Gamma(y+1)}{x y \Gamma(x+y+2)} \ dx \ dy$$ 
that is obtained by manipulating 
$$\int_0^1 \operatorname{li}(x) \operatorname{li}(1-x) \ dx$$
 Apparently nothing seems to work, but I bet you can do much more than me. So, what
tools would you recommend me to try?
 A: To compute
$$\int_{0}^{1}\operatorname{li}(x)\operatorname{li}(1-x)\,dx$$
i'd rather exploit the fact that the logarithmic integral function admits a representation in terms of the shifted Legendre polynomials, since:
$$\begin{eqnarray*}\int_{0}^{1}\operatorname{li}(x)\,L_n(2x-1)\,dx &=& \frac{1}{n!}\int_{0}^{1}\operatorname{li}(x)\frac{d^n}{dx^n}(x^2-x)^n\,dx\\&=&\frac{1}{n!}\left.\frac{L_{n+1}(2x-1)-L_{n-1}(2x-1)}{4n+2}\operatorname{li}(x)\right|_{0}^{1}\\&-&\frac{1}{(4n+2)n!}\int_{0}^{1}\frac{L_{n+1}(2x-1)-L_{n-1}(2x-1)}{\log x}\,dx.\end{eqnarray*}$$
So, if we set $L_n(2x-1)=Q_n(x)$ for shortness, we have:
$$\operatorname{li}(x) = -\log 2-\sum_{n=1}^{+\infty} c_n Q_n(x), $$
where:
$$ c_n = \frac{1}{2n!}\int_{0}^{1}\frac{Q_{n+1}(x)-Q_{n-1}(x)}{\log x}\,dx \tag{1}$$
is a linear combination of logarithms of integers by Frullani's theorem and:
$$\operatorname{li}(1-x) = -\log 2-\sum_{n=1}^{+\infty}(-1)^n c_n Q_n(x), $$
giving:
$$\int_{0}^{1}\operatorname{li}(x)\operatorname{li}(1-x)\,dx=\log^2 2+\sum_{n=1}^{+\infty}\frac{(-1)^n c_n^2}{2n+1}.\tag{2}$$
