Step 1: If $x_n = \tan(\phi)$ for some $\phi $ then
$$
x_{n+1} = -\frac{2 \tan(\phi)}{1- \tan^2(\phi)} = - \tan(2 \phi) = \tan(-2\phi) \, .
$$
Therefore, if we define
$$
\alpha = \arctan(x_1), \, \beta = \arctan(y_1), \, \gamma = \arctan(z_1) \, ,
$$
then
$$
x_n = \tan((-2)^{n-1} \alpha), \, y_n = \tan((-2)^{n-1} \beta), \, z_n = \tan((-2)^{n-1} \gamma)
$$
for all $n$, as long as the sequences are defined, i.e. as long as no denominator in the sequence definitions becomes zero.
Step 2: Generally,
$$
\tan (u) + \tan (v) + \tan (w) = \tan (u) \tan(v) \tan(w)
$$
holds if and only if $u+v+w$ is an integral multiple of $\pi$, see for example Prove the identity if x+y+z=xyz on AoPS.
Step 3: The initial values $x_1 = 2$, $y_1 = 4$, $z_1 = 6/7$ satisfy $x_1 + y_1 + z_1 = x_1 y_1 z_1$. It follows that $\alpha + \beta + \gamma$ is an integral multiple of $\pi$. But then also
$$
(-2)^{n-1} \alpha + (-2)^{n-1} \beta + (-2)^{n-1} \gamma = k \pi \text{ for some }k \in \Bbb Z
$$
which in turn implies that
$$
x_n + y_n + z_n = x_n y_n z_n
$$
as long as the sequences are defined. It follows directly from the sequence definitions that the right-hand side is never zero, therefore $x_n + y_n + z_n$ can not be zero either.
Remark: One can also verify directly from the recursive definition that
$$
x_{n+1} + y_{n+1} + z_{n+1} - x_{n+1} y_{n+1} z_{n+1} =
\frac{2(1-x_n y_n - x_n z_n - y_n z_n)(x_n + y_n +z_n -x_ny_n z_n)}{(x_n^2-1)(y_n^2-1)(z_n^2-1)}
$$
so that
$$
x_n + y_n +z_n = x_ny_n z_n \implies x_{n+1} + y_{n+1} + z_{n+1} = x_{n+1} y_{n+1} z_{n+1} \, .
$$
I suspected something deeper was going on. Indeed, if we just substitute $L$ into the recursions, we get trivial equations. So clearly we can't work our way out that easily!
I've often seen situations where invariants come in useful. Given a sequence $x_n$, a function $f$ such that $f(x_n) = f(x_{n+1})$ for all $n$ is called an invariant of the sequence.
The key thing, is that if $f$ is continuous then $f(x_0) = f(x_1)=\ldots = f(\lim_{n \to \infty} x_n)$, and therefore we may be able to use the nature of preimages of $f$ to calculate our answer.
The first thing I did was write Python code to simulate the situation. I just took some arbitrary $a,b,c$ here.
from math import *
def checkone() :
a = 3
b = 7
c = sqrt(61)
for i in range(100) :
d = a
e = b
f = c
a = sqrt(d*(f+e-d))
b = sqrt(e*(f+d-e))
c = sqrt(f*(d+e-f))
print(str([a,b,c]) + ","+str(a**4+b**4+c**4))
checkone()
where I used the $a^4+b^4+c^4$ term to predict that something was going on with the fourth powers of the sequence. I noticed this fourth-power pattern when I saw that we would need to play with quadratics of $x_n^2$, which due to the nature of the RHS (the fact that $x_{n+1}^2 = sth. x^2 + sth. x_ny_n $ etc.) meant that I'd have to go a square further.
Call $(L(a,b,c),L(a,b,c),L(a,b,c))$ as the limit obtained from the above. Following some playing around, I first came to the conclusion that :
$$
3L^4(a,b,\sqrt{a^2+b^2}) = 4a^2b^2
$$
(note : I actually used the homogeniety of the relation and realized that it was enough to take $a=1$, so this simplified things) and then by perturbation to the conclusion:
$$
3L^4(a,b,\sqrt{a^2+b^2-c}) = 4a^2b^2-c^2
$$
This then tells me what the limit should be , following a change of subject:
$$
3L^4(a,b,x) = 2(a^2b^2+b^2x^2+x^2a^2) - a^4-b^4-x^4 \\
$$
Now comes the wonderful idea : That quantity, right there on the RHS, is actually invariant.
Suppose that $x_n,y_n,z_n$ are defined as given. Then :
$$
2(x_{n+1}^2y_{n+1}^2+x_{n+1}^2z_{n+1}^2+y_{n+1}^2z_{n+1}^2) - x_{n+1}^4-y_{n+1}^4-z_{n+1}^4\\ = 2(x_n^2y_n^2+y_n^2z_n^2+z_n^2x_n^2)-x_n^4-y_n^4-z_n^4
$$
This is a non-routine but fairly ok check to make. So you substitute the recursion where required and make sure things simplify.
Once we make this check, we note that we can, obviously by induction, obtain the statement :
$$
2(x_{n}^2y_{n}^2+x_{n}^2z_{n}^2+y_{n}^2z_{n}^2) - x_{n}^4-y_{n}^4-z_{n}^4\\ = 2(x_0^2y_0^2+y_0^2z_0^2+z_0^2x_0^2)-x_0^4-y_0^4-z_0^4
$$
for all $n \geq 1$. Now, the splendid work of the OP has proved that the sequence converges to a limit $(\alpha,\alpha,\alpha)$.
But now, we may take the limit on both sides of the recursion. Indeed, note that if we do so then we get :
$$
3L(x_0,y_0,z_0)^4 = 2(x_0^2y_0^2+y_0^2z_0^2+z_0^2x_0^2)-x_0^4-y_0^4-z_0^4
$$
and therefore :
$$
\boxed{
L(x_0,y_0,z_0) = \sqrt[4]{\frac{2(x_0^2y_0^2+y_0^2z_0^2+z_0^2x_0^2)-x_0^4-y_0^4-z_0^4}{3}}
}
$$
which gives the correct answer for the case above.
So find invariants, is the lesson. Don't ask me how I'd do this without computing power, though!
Best Answer
The fractional part is not dense when $b$ is a square. Indeed, if $b = j^2$, then
$$x_{n} = 2\sqrt{j^{2n}(j^{2n+2}-1)},$$ which is very slightly smaller than the integer $2j^{2n+1}$. In fact,
$$ \begin{aligned} 0 < 2j^{2n+1} -2\sqrt{j^{2n}(j^{2n+2}-1)} &= \frac{(2j^{2n+1} -2\sqrt{j^{2n}(j^{2n+2}-1)})(2j^{2n+1} +2\sqrt{j^{2n}(j^{2n+2}-1)})}{2j^{2n+1} +2\sqrt{j^{2n}(j^{2n+2}-1)}}\\ &= \frac{4j^{2n}}{2j^{2n+1} +2\sqrt{j^{2n}(j^{2n+2}-1)}}\\ &=\frac{2}{j + \sqrt{j^2-j^{-2n}}}\\ &< 1. \end{aligned} $$
Thus, if $b$ is a square, then $\{x_{n}\}$ is strictly increasing, and furthermore $\{x_{n}\} \rightarrow 1 - 1/j$.