A triple integral dancing in the unit cube Straight integration seems pretty tedious and difficult, and I guess the symmetry may open some new ways of which I'm not aware. What would your idea be?
Calculate
$$\int_0^1 \int_0^1 \int_0^1 \frac{x^2}{\sqrt{x^2+1} \left(x^2-y^2\right) \left(x^2-z^2\right)}+\frac{y^2}{\sqrt{y^2+1} \left(y^2-x^2\right) \left(y^2-z^2\right)}+\frac{z^2}{\sqrt{z^2+1} \left(z^2-x^2\right) \left(z^2-y^2\right)} \, dx \ dy \ dz.$$
A 300 points bounty moment: After 2 years and 10 months since the problem was posed no full solution was provided yet. Is it possible to find a slick solution (like a bolt of lightning)? Good luck!
 A: A note: The integral is only apparently improper. Substituting
$$x={1\over2}\left(u-{1\over u}\right)\quad(1\leq u\leq1+\sqrt{2}),\quad\sqrt{1+x^2}={1\over2}\left(u+{1\over u}\right),\quad dx={1\over2}\left(1+{1\over u^2}\right)\>du\ ,$$
and similarly for the other variables leads to an integrand with a polynomial in the numerator and
$$(u+v)(1+uv)(v+w)(1+vw)(w+u)(1+wu)$$
in the denominator. Mathematica computes the following numerical value of the integral (in $(u,v,w)$-terms as well as in the original form):
$$-0.287162469\ .$$
A: Let us set $A=x^2,B=y^2,C=z^2$ and:
$$ S_n(A,B,C) = \frac{A^n}{(A-B)(A-C)}+\frac{B^n}{(B-A)(B-C)}+\frac{C^n}{(C-A)(C-B)} .\tag{1}$$
By partial fraction decomposition, it is straightforward to check that $S_0=S_1=0$ and $S_2=1$.
For every $n>2$ induction gives:
$$ S_n(A,B,C) = h_{n-2}(A,B,C)\tag{2} $$
where $h_k$ is a complete homogeneous symmetric polynomial, satisfying:
$$ \sum_{k\geq 0}h_k(A,B,C)t^k=\frac{1}{(1-A t)(1-B t)(1- Ct)}.\tag{3} $$
Assume now that the Taylor series of $g(u)=\frac{u}{\sqrt{1+u}}$ around $u=0$ is given by $\sum_{n\geq 0}g_n u^n$.
Our integral is given by:
$$ \iiint_{(0,1)^3}\sum_{n\geq 0}g_n\,S_n(x^2,y^2,z^2)\,d\mu=g_2+\iiint_{(0,1)^3}\sum_{n\geq 3}g_n\cdot h_{n-2}(x^2,y^2,z^2)\,d\mu\tag{4} $$
Now there are two possible approaches. The first approach is to look for  a linear operator that maps $t^k$ to $g_k$, then exploit $(3)$ in order to write the RHS of $(4)$ as a simpler integral. The second approach is to apply a double-counting argument on the monomials appearing in the RHS of $(4)$. Since $\iiint_{(0,1)^3}x^{2a}y^{2b}z^{2c}\,d\mu=\frac{1}{(2a+1)(2b+1)(2c+1)}$, we get:
$$ \iiint_{(0,1)^3}\sum_{n\geq 3}g_n\cdot h_{n-2}(x^2,y^2,z^2)\,d\mu =\\= \sum_{\substack{a,b,c\geq 0\\(a,b,c)\neq (0,0,0)}}\frac{g_{a+b+c+2}}{(2a+1)(2b+1)(2c+1)}\tag{5} $$
and the last series can be further simplified by computing:
$$ l_n=\sum_{\substack{a,b,c\geq 0\\a+b+c=n}}\frac{1}{(2a+1)(2b+1)(2c+1)}\tag{6} $$
then computing $\sum_{n\geq 1}l_n g_n$, that is obviously related with the integral: $$\iiint_{(0,1)^3}g(x^2 y^2 z^2)\,d\mu = \iiint_{(0,1)^3}\frac{x^2 y^2 z^2}{\sqrt{1+x^2 y^2 z^2}}\,d\mu.\tag{7}$$

2018 update. Following the first approach, the given integral equals the opposite of
$$\iiint_{(0,1)^3}\frac{2}{\pi}\int_{0}^{\pi/2}\frac{\sin^2(\theta)\,d\theta}{(1+x^2\sin^2\theta)(1+y^2\sin^2\theta)(1+z^2\sin^2\theta)}\,dx\,dy\,dz$$
which by Fubini's theorem boils down to
$$ \frac{2}{\pi}\int_{0}^{\pi/2}\frac{\arctan^3(\sin\theta)}{\sin\theta}\,d\theta = \frac{2}{\pi}\int_{0}^{\pi/2}\frac{\arctan^3(\cos\theta)}{\cos\theta}\,d\theta=\frac{4}{\pi}\int_{0}^{1}\frac{\arctan^3\left(\frac{1-t^2}{1+t^2}\right)}{1-t^2}\,dt $$
then to
$\frac{2}{\pi}\int_{0}^{1}\frac{\arctan^3\left(\frac{1-t}{1+t}\right)}{\sqrt{t}(1-t)}\,dt = \frac{2}{\pi}\int_{0}^{1}\frac{\left(\frac{\pi}{4}-\arctan t\right)^3}{\sqrt{t}(1-t)}\,dt = \frac{2}{\pi}\int_{0}^{+\infty}\left(\frac{\pi}{4}-\arctan\tanh^2\frac{u}{2}\right)^3\,du$ which is related to the Gudermannian function. Probably it can be expressed in a closed form through high order polylogarithms / Euler sums with high weight. Accurate numerical approximation are pretty simple since the last integrand function essentially behaves like $\frac{\pi^3}{64}\exp\left(-\frac{3}{\pi}u^2\right)$ on $\mathbb{R}^+$.
We may also notice that the two approaches are pretty obviously equivalent, since $l_n$ appearing in $(6)$ is clearly given by a coefficient of the MacLaurin series of $\text{arctanh}^3(z)$.
