Solution (source code)

= Solution
{c}

Choose the first <orthonormal basis> <vector> to be $\mathbf e_1=\mathbf a$. If $\mathbf b$ is not parallel to $\mathbf a$, choose $\mathbf e_2$ along $\mathbf b-(\mathbf a\cdot\mathbf b)\mathbf a$; otherwise choose any <unit vector> <orthogonal> to $\mathbf a$. Complete to an <orthonormal basis> with $\mathbf e_3$. The <matrix> with rows $\mathbf e_i^T$ is an <orthogonal matrix>, and expressing the three <unit vectors> in this <basis> gives
$$
\mathbf a'=(1,0,0)^T,\qquad \mathbf b'=(x,y,0)^T,\qquad
\mathbf c'=(z,w\cos\theta,w\sin\theta)^T.
$$
The last two components of $\mathbf c'$ have <polar coordinates> of the displayed form; we may choose $w\geq0$. If they vanish, $\theta$ is arbitrary.

An <orthogonal transformation> preserves <dot products>. Taking <dot products> with $\mathbf a'$, and then using the unit <Euclidean norms>, gives
$$
\boxed{x=\cos\gamma,\qquad z=\cos\beta,\qquad |y|=\sin\gamma,\qquad |w|=\sin\beta.}
$$
The remaining <dot product> gives
$$
\cos\alpha=xz+yw\cos\theta.
$$
An <orthogonal matrix> has <determinant> $\pm1$, so it preserves the absolute <scalar triple product>. In the new <coordinates> this <scalar triple product> is $yw\sin\theta$. Consequently
$$
\begin{aligned}
| (\mathbf a\wedge\mathbf b)\cdot\mathbf c|^2
&=y^2w^2-(yw\cos\theta)^2\\
&=(1-x^2)(1-z^2)-(\cos\alpha-xz)^2\\
&=1-\cos^2\alpha-\cos^2\beta-\cos^2\gamma
+2\cos\alpha\cos\beta\cos\gamma.
\end{aligned}
$$
Taking the nonnegative <square root> proves
$$
\boxed{|(\mathbf a\wedge\mathbf b)\cdot\mathbf c|
=\sqrt{1-(\mathbf a\cdot\mathbf b)^2-(\mathbf b\cdot\mathbf c)^2-(\mathbf c\cdot\mathbf a)^2+2(\mathbf a\cdot\mathbf b)(\mathbf b\cdot\mathbf c)(\mathbf c\cdot\mathbf a)}.}
$$
No division by $y$ or $w$ was used, so the proof includes parallel <vectors>. Equivalently, the expression under the <square root> is the <Gram determinant> of the three <unit vectors>.