Solution (source code)

= Solution

Use <metric signature> $(+---)$ and <four-velocities> $U_i^a=\gamma_i c(1,\boldsymbol\beta_i)$, with $\gamma_i=(1-\beta_i^2)^{-1/2}$. Their <Lorentz invariant> <inner product> is
$$
\frac{U_1\cdot U_2}{c^2}
=\gamma_1\gamma_2(1-\boldsymbol\beta_1\cdot\boldsymbol\beta_2).
$$
In the <rest frame> of either particle this same <scalar> is the relative <Lorentz factor> $\gamma_{\rm rel}$. Hence
$$
\beta_{\rm rel}^2
=1-\frac{(1-\beta_1^2)(1-\beta_2^2)}
{(1-\boldsymbol\beta_1\cdot\boldsymbol\beta_2)^2}.
$$
Let $d=\boldsymbol\beta_1\cdot\boldsymbol\beta_2$. The numerator after combining terms is
$$
(1-d)^2-(1-\beta_1^2)(1-\beta_2^2)
=\beta_1^2+\beta_2^2-2d-\beta_1^2\beta_2^2+d^2.
$$
Use $|\boldsymbol\beta_1\times\boldsymbol\beta_2|^2=\beta_1^2\beta_2^2-d^2$ to obtain
$$
\boxed{\beta_{\rm rel}^2=
\frac{|\boldsymbol\beta_1-\boldsymbol\beta_2|^2-
|\boldsymbol\beta_1\times\boldsymbol\beta_2|^2}
{(1-\boldsymbol\beta_1\cdot\boldsymbol\beta_2)^2}.}
$$
The two <rest frames> agree on the relative speed, although their <velocity vectors> are expressed in different spatial frames. Checks include equal velocities giving zero, collinear velocities giving $|\beta_1-\beta_2|/|1-\beta_1\beta_2|$, and low speeds reducing to the Galilean relative speed. For timelike particles the <four-velocity> <inner product> also ensures $0\le\beta_{\rm rel}<1$.