<s>[INST] How to solve this 2-D integration? How to solve 
$$I=\int_{[0,1]^2}\frac{dxdy}{(1+xy)(1+x^2)}$$
I've tried using the diffeomorphism $(x,y)=(u,v/u)$ from $\text{int}\{(u,v)\mid 0\le u\le 1,0\le v\le u\}$ to $\text{int}[0,1]^2$, but since the Jaccobian determinant becomes unbounded at the boundary, the transformation leads to an improper integral, which I think makes it invalid. 
If I integrate directly over the original domain, the result will be too complicated to calculate. 
I can't come up with another diffeomorphism. Is there another method to go on?
 [/INST] I don't know if this is the slickest way, but at least it works. First decompose into partial fractions with respect to $x$:
$$
\frac{1}{(1+xy)(1+x^2)} = \frac{1}{1+y^2} \left( \frac{y^2}{1+xy} + \frac{1-xy}{1+x^2} \right)
.
$$
Rewrite this (using $y^2=(1+y^2)-1$ as)
$$
\frac{1}{1+xy}
- \frac{1}{(1+xy)(1+y^2)}
+ \frac{1}{1+x^2} \cdot \frac{1}{1+y^2}
- \frac{x}{1+x^2} \cdot \frac{y}{1+y^2}
.
$$
The last two terms are easy to integrate over the square, since $x$ and $y$ separate.
The integral of the second term equals the integral that we're trying to compute (call it $I$), by symmetry.
The integral of the first term is trickier, but we can expand $(1+xy)^{-1}=1 - x y  + x^2 y^2 - x^3 y^3 + \dotsb$, which upon termwise integration gives the known series $1-1/2^2+1/3^2-1/4^2+\dotsb = \pi^2/12$. One could probably also get this in the same way as the integral of $(1-xy)^{-1}$ over the unit square is computed in Chapter 8 of Proofs from THE BOOK by Aigner & Ziegler, namely by setting $x=u-v$ and $y=u+v$.
All combined, we get
$$
I = \frac{\pi^2}{12} - I + \left( \frac{\pi}{4} \right)^2 - \left( \frac{\ln 2}{2} \right)^2
,
$$
so
$$
I=\frac{7 \pi^2}{96} - \frac{\ln^2 2}{8}
.
$$
</s>