Evaluate $ \int_0^{\infty} \int_0^{\infty} \frac{1}{(1+x)(1+y)(x+y)}\,dy dx$ 
Evaluate $ \int_0^{\infty} \int_0^{\infty} \frac{1}{(1+x)(1+y)(x+y)}\,dy dx$

We can calculate $\int \frac{1}{(1+x)(x+y)} \,dx= \frac{\ln|x+1| - \ln|x+y|}{y-1}$ but that doesn't make it quite easier. Is there a trick to compute this easier?
 A: CW because some ideas were taken from other posts, which I have linked to.
I'll integrate first in $y$ and then reduce the problem to a known integral. Using partial fractions, we have$$
\frac{1}{(1+y)(x+y)} = \frac{1}{(x-1)}\left(\frac{1}{1+y}-\frac{1}{x+y}\right)
$$Thus
$$
\int_0^{\infty}\int_0^{\infty}\frac{1}{(1+x)(1+y)(x+y)} \,dydx 
$$
$$
=\int_0^{\infty}\frac{1}{x^2-1}\int_0^{\infty}\frac{1}{1+y}-\frac{1}{x+y} \,dydx 
$$
$$
=\int_0^{\infty}\frac{1}{x^2-1}\left(\left.\log\left|\frac{1+y}{x+y}\right|\right|_0^{\infty}\right)\,dx 
$$
$$
=\int_0^{\infty}\frac{1}{x^2-1}\left(-\log\left|\frac{1}{x}\right|\right)\,dx 
$$
$$
=\int_0^{\infty}\frac{\log|x|}{x^2-1}\,dx 
$$Break into $[0,1)$ and $(1,\infty)$ (the function has a removable singularity at $x=1$, is integrable near zero by the MVT and at infinity by direct comparison, so we're good) and enforce the substitution $z=1/x$. You'll see that they are equal; hence it suffices to integrate over $(0,1)$ and double it. This allows us to use the Maclaurin series of log, convergent on $(0,1)$:
$$
\int_0^{\infty}\frac{\log|x|}{x^2-1}\,dx =2 \int_0^1 \frac{\log(x)}{(x+1)(x-1)}\,dx
$$It's a nice exercise to work out that this last integral gives $\pi^2/8$ by reducing it to $\sum_{n\ge 1}(2n-1)^{-2}$. More generally, we have $\int _0^{\infty} \frac{\log(x)}{x^2+\alpha^2}\,dx = \frac{\pi\log \alpha}{2\alpha}$, and interpreting $\log(i) = \pi i/2$ gives us $\pi^2/4$, as above.
A: $\newcommand{\bbx}[1]{\,\bbox[15px,border:1px groove navy]{\displaystyle{#1}}\,}
 \newcommand{\braces}[1]{\left\lbrace\,{#1}\,\right\rbrace}
 \newcommand{\bracks}[1]{\left\lbrack\,{#1}\,\right\rbrack}
 \newcommand{\dd}{\mathrm{d}}
 \newcommand{\ds}[1]{\displaystyle{#1}}
 \newcommand{\expo}[1]{\,\mathrm{e}^{#1}\,}
 \newcommand{\ic}{\mathrm{i}}
 \newcommand{\mc}[1]{\mathcal{#1}}
 \newcommand{\mrm}[1]{\mathrm{#1}}
 \newcommand{\on}[1]{\operatorname{#1}}
 \newcommand{\pars}[1]{\left(\,{#1}\,\right)}
 \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}}
 \newcommand{\root}[2][]{\,\sqrt[#1]{\,{#2}\,}\,}
 \newcommand{\totald}[3][]{\frac{\mathrm{d}^{#1} #2}{\mathrm{d} #3^{#1}}}
 \newcommand{\verts}[1]{\left\vert\,{#1}\,\right\vert}$
\begin{align}
{\cal J} &\equiv
\bbox[5px,#ffd]{\int_{0}^{\infty}\int_{0}^{\infty}
{\dd x\,\dd y \over
\pars{1 + x}\pars{1 + y}\pars{x + y}}}:\ {\Large ?}.
\end{align}
With the variable changes $\ds{x \equiv 1/a - 1}$ and
$\ds{y \equiv 1/b - 1}$ the above integral becomes
\begin{align}
{\cal J} & = \int_{0}^{1}\int_{0}^{1}
{\dd a\,\dd b \over a + b - 2ab} =
\int_{0}^{1}\ln\pars{1 - b \over b}
{\dd b \over 1 - 2b}
\end{align}
With $\ds{\pars{~1 - 2b = t \implies
b = {1 - t \over 2}~}}$:
\begin{align}
{\cal J} & = \underbrace{\int_{0}^{-1}{\ln\pars{1 - t} \over t}\,\dd t}_{\ds{\pi^{2} \over 12}}\ -\
\underbrace{\int_{0}^{1}{\ln\pars{1 - t} \over t}\,\dd t}_{\ds{-\,{\pi^{2} \over 6}}}
\\[5mm] & = \bbx{\large{\pi^{2} \over 4}} \\&&
\end{align}
