Solution (source code)

= Solution

Contract the <Levi-Civita symbol> identity $\varepsilon_{ijk}\varepsilon_{i\ell m}=\delta_{j\ell}\delta_{km}-\delta_{jm}\delta_{k\ell}$. Using <Einstein summation convention>,
$$
\begin{aligned}
|\mathbf a\times\mathbf b|^2
&=\varepsilon_{ijk}a_jb_k\varepsilon_{i\ell m}a_\ell b_m\\
&=(\delta_{j\ell}\delta_{km}-\delta_{jm}\delta_{k\ell})a_jb_ka_\ell b_m\\
&=|\mathbf a|^2|\mathbf b|^2-(\mathbf a\cdot\mathbf b)^2.
\end{aligned}
$$
This <Lagrange identity for the cross product> includes zero and parallel vectors without a separate division by their norms.