Let $A := \{ (x, y, z) \in \mathbb{R}^{3} : 0 < x < 1, 0 < y < 1, 0 < z \} $.
Show that the function $f: A \to \mathbb{R}$ defined by $$f(x, y, z) := \frac{1}{(1+x^2z^2)(1+y^2z^2)}$$ is integrable in $A$ and calculate the integral $\int_{A} f$.
My attempt: Showing the integral exists is easy using the Comparison test for improper integrals. Now, for computing the integral first I've used Fubini's theorem, so I got that $$\int_{A} f = \int_{0}^{\infty} \left( \frac{\arctan(z)}{z} \right)^{2} \: dz \,,$$ which is not an elementary integral. Then I thought about using that $f$ is an even function, and I got that $$\int_{A} f = \frac{1}{8} \int_{D} f\,,$$ where $D:= \{ (x, y, z): -1 < x <1, -1 < y <1 \}$ and tried to use a change of variables to cylindrical coordinates...what doesn't make things easier. I also thought about non-linear changes of variables such as $u := xz$, $v := yz$, $z = z$...
But I am no longer able to solve the integral.
It is more convenient to calculate the integral with respect to $z$ first and then with respect to $x$ and $y$.
We have $$\int_{0}^{\infty}\frac{1}{(1+x^2 z^2)(1+y^2 z^2)} dz=\frac{x \operatorname{arctg}xz-y\operatorname{arctg}yz}{x^2-y^2}\Big|_0^\infty=\frac{\pi}{2(x+y)}$$ Now we find $$\int_0^{1}\frac{\pi}{2(x+y)}dy=\frac{\pi}{2}(\ln(1+x)-\ln(x))$$ and finally $$\int_{0}^{1}\frac{\pi}{2}(\ln(1+x)-\ln(x))dx=\frac{\pi}{2}\left((1+x)\ln(\frac{1+x}{x})-x \ln(x) \right)\Big|_0^1=\pi \ln(2)$$
since $\int \ln(x) dx=x(\ln(x)-1)$. $\;$ Approximately $\pi \ln(2)\approx 2.17759$.
Alternatively a few tricks how to calculate $\;\int_0^\infty \frac{\operatorname{arctg}^2 z}{z^2} dz \;$ one can find here Integral $\int_0^{\infty}\frac{\operatorname{arctg}^2x}{x^2}dx$.