Solution (source code)

= Solution

For $v=|\mathbf v|>0$, decompose the <velocity> as $\mathbf u=\mathbf u_\parallel+\mathbf u_\perp$, parallel and perpendicular to $\mathbf v$. A rotation-free <Lorentz boost> gives the <velocity-addition formula>
$$
\mathbf w_\parallel=\frac{\mathbf u_\parallel-\mathbf v}{D},\qquad
\mathbf w_\perp=\frac{\mathbf u_\perp}{\gamma_vD},\qquad
D=1-\mathbf v\cdot\mathbf u/c^2.
$$
These two components are orthogonal, so
$$
\begin{aligned}
D^2|\mathbf w|^2
&=|\mathbf u_\parallel-\mathbf v|^2+(1-v^2/c^2)|\mathbf u_\perp|^2\\
&=|\mathbf u-\mathbf v|^2-\frac{v^2|\mathbf u_\perp|^2}{c^2}\\
&=|\mathbf u-\mathbf v|^2-\frac{|\mathbf v\times\mathbf u|^2}{c^2}.
\end{aligned}
$$
The last equality is <Lagrange's identity> for the <cross product>. For physical <subluminal speeds>, $D>0$, giving
$$
\boxed{|\mathbf w|=\frac{\sqrt{(\mathbf v-\mathbf u)\cdot(\mathbf v-\mathbf u)-|\mathbf v\times\mathbf u|^2/c^2}}{1-\mathbf v\cdot\mathbf u/c^2}.}
$$
The $\mathbf v=0$ case follows directly or by continuity. An independent form useful for checks is $1-|\mathbf w|^2/c^2=(1-u^2/c^2)(1-v^2/c^2)/D^2$, obtained from the transformed <four-velocity> normalization.