numerical calculation of 3d integral $\iiint_{[-1,1]^3} \frac{dx dy dz}{r^2} $ I have to calculate the integral
$$I=\iiint_{[-1,1]^3} \frac{dx dy dz}{r^2} $$
where $r^2 = x^2 + y^2 + z^2 $.
My strategy is to split the integral domain into two parts: the unit sphere and the supplement domain. From the unit sphere, I get analytically $4\pi $; for the supplement, I used Monte Carlo.  The problem is that the Monte Carlo part is not very efficient, while I want to get the integral to 8 digits.
Could anyone suggest a better scheme?
 A: Per symmetry, the integral over $[-1,1]^3$ is $8$ times that over $[0,1]^3$; the latter is $6$ times the integral over the tetrahedron $0<z<y<x<1$. Thus
$$\iiint_{[-1,1]^3} \frac{dx dy dz}{r^2}= 48J$$
where
\begin{align}
J&=\int_0^1 \int_0^x \int_0^y \frac{1}{x^2+ y^2+z^2}\ \overset{z=yt}{dz}\ dy \ dx\\
&= \int_0^1 \int_0^1\int_0^x \frac{y}{x^2+ y^2(1+t^2)}\ \overset{y=x \sinh v}{dy\ \ dx } \ dt\\
 &= \int_0^{\sinh^{-1}1}\int_0^1  \frac{\tanh v }{1+ \tanh^2 v \ t^2} dt\ dv
= \int_0^{\sinh^{-1}1} \tan^{-1}(\tanh v)
\overset{u=e^{-2v}}{dv}\\
&=\ \frac12 \int_{(\sqrt2-1)^2}^1 \frac1u{\tan^{-1}\frac{1-u}{1+u}}\ {du}
=\frac12 \int_{(\sqrt2-1)^2}^1 \frac{\frac\pi4-\tan^{-1}u}{u}{du}\\
&=-\frac\pi8\ln(\sqrt2-1)^2-\frac12\int_0^1 \frac{\tan^{-1}u}udu
 + \frac12\int_0^{(\sqrt2-1)^2} \frac{\tan^{-1}u}udu \\
&=- \frac\pi4\ln(\sqrt2-1)-\frac12G+\frac12 \Im\text{Li}_2(i(\sqrt2-1)^2)
\end{align}
A: $$\iiint_{[-1,1]^3}\frac{d\mu}{x^2+y^2+z^2} = \iint_{[-1,1]^2}\frac{2}{\sqrt{x^2+y^2}}\arctan\left(\frac{1}{\sqrt{x^2+y^2}}\right)\,dx\,dy $$
equals
$$ 8\iint_{[0,1]^2}\frac{1}{\sqrt{x^2+y^2}}\arctan\left(\frac{1}{\sqrt{x^2+y^2}}\right)\,dx\,dy $$
or
$$ 16\iint_{0\leq y\leq x\leq 1}\frac{1}{\sqrt{x^2+y^2}}\arctan\left(\frac{1}{\sqrt{x^2+y^2}}\right)\,dx\,dy \stackrel{y\mapsto tx}{=} 16\iint_{[0,1]^2}\frac{1}{\sqrt{1+ t^2}}\arctan\left(\frac{1}{x\sqrt{1+ t^2}}\right)\,dx\,dt $$
or
$$ 16\int_{0}^{1}\left(\frac{\pi}{2\sqrt{1+t^2}}-\frac{\arctan\sqrt{1+t^2}}{\sqrt{1+t^2}}+\frac{\log(2+t^2)}{2(1+t^2)}\right)\,dt $$
or
$$ 8\pi\log(1+\sqrt{2})+16\int_{0}^{1}\left(-\frac{\arctan\sqrt{1+t^2}}{\sqrt{1+t^2}}+\frac{\log(2+t^2)}{2(1+t^2)}\right)\,dt. $$
The last integral can be written in terms of polylogarithms, or just numerically evaluated through Newton-Cotes or Gaussian quadrature. The outcome, of course, is larger than $4\pi$ and it equals $\color{blue}{15.348248444887\ldots}$
