Solution (source code)

= Solution

The nonparallel direction <unit vectors> have nonzero <cross product>. Let
$$
\mathbf n=\frac{\widehat{\mathbf b}\times\widehat{\mathbf d}}
{|\widehat{\mathbf b}\times\widehat{\mathbf d}|}.
$$
For any joining <vector> $\mathbf s=\mathbf c-\mathbf a+\mu\widehat{\mathbf d}-\lambda\widehat{\mathbf b}$, its normal component is fixed:
$$
\mathbf s\cdot\mathbf n=(\mathbf c-\mathbf a)\cdot\mathbf n.
$$
Thus its length is at least the absolute value of that component. Since the two line directions span $\mathbf n^\perp$, suitable $\lambda,\mu$ cancel the entire component in that <plane>, attaining the bound. Therefore the <distance between skew lines>, including the zero-distance intersecting case, is
$$
\boxed{d_{\min}=\frac{|(\mathbf c-\mathbf a)\cdot
(\widehat{\mathbf b}\times\widehat{\mathbf d})|}
{|\widehat{\mathbf b}\times\widehat{\mathbf d}|}.}
$$
The numerator is a <scalar triple product>. The nonparallel hypothesis is essential to the denominator and to the spanning argument.