compute the integral $\int_0^1 \int_0^1 \int_0^1 \frac{1}{(1+x^2+y^2+z^2)^2} dxdydz$ 
Determine, with justification, the value of the integral $\int_0^1 \int_0^1 \int_0^1 \frac{1}{(1+x^2+y^2+z^2)^2} dxdydz$.

I tried converting this integral to cylindrical coordinates with $r = \sqrt{x^2 + y^2}$ ranging from $0$ to $\sqrt{2}$, $0\leq \theta \leq \pi/2, 0\leq z \leq 1,$ where $\theta $ is such that $x= r\cos\theta, y = r\sin\theta.$ However, this seems to lead to an incorrect result. Which bounds have I gotten wrong? Also, it seems that the integral over the unit cube equals twice the integral over the region defined by $0\leq z\leq 1, 0\leq x\leq 1, 0\leq y\leq x,$ but I'm not sure why. The result should be $\frac{\pi^2}{32},$ which is basically what WolframAlpha outputs.
Using spherical coordinates seems to make the integration more complicated due to the integration factor.
 A: Using the identity ($q > 0$)
$$\frac{1}{q^2} = \int_0^\infty t\mathrm{e}^{-qt} \,\mathrm{d} t,$$
we have
\begin{align*}
 I &= \int_0^1\int_0^1 \int_0^1
 \int_0^\infty t\mathrm{e}^{-(1 + x^2 + y^2 + z^2)t} \,\mathrm{d} t \,\mathrm{d}x \mathrm{d}y \mathrm{d}z\\
 &= \int_0^\infty t\mathrm{e}^{-t} \left(\int_0^1 \mathrm{e}^{- x^2t}\mathrm{d} x\right)^3\mathrm{d}t\\
 &= \int_0^\infty t\mathrm{e}^{-t} \left(\frac{1}{\sqrt{t}}\int_0^{\sqrt{t}} \mathrm{e}^{- u^2}\mathrm{d} u\right)^3\mathrm{d}t\\
 &= \int_0^\infty \frac{1}{\sqrt{t}}\mathrm{e}^{-t} \left(\int_0^{\sqrt{t}} \mathrm{e}^{- u^2}\mathrm{d} u\right)^3\mathrm{d}t \\
 &= \frac12\left(\int_0^{\sqrt{t}} \mathrm{e}^{- u^2}\mathrm{d} u\right)^4\Big\vert_0^\infty \\
 &= \frac12 \left(\int_0^{\infty} \mathrm{e}^{- u^2}\mathrm{d} u\right)^4\\
 &= \frac12\left(\frac12\sqrt{\pi}\right)^4\\
 &= \frac{\pi^2}{32}
\end{align*}
where we have used $\int_{-\infty}^\infty \mathrm{e}^{-u^2}\mathrm{d} u = \sqrt{\pi}$,
and $(\int_0^{\sqrt{t}} \mathrm{e}^{- u^2}\mathrm{d} u)'
= \frac{1}{2\sqrt{t}}\mathrm{e}^{-t}$.
