It's a one-liner in quaternion algebra: vu-|vu|

followed by normalization. Proof can be carried out completely in quaternion algebra.

An interesting degenerate case is when u=-v. In that case, another expression can be used: PERP(v)(vu+|vu|) {where PERP(v) means an arbitrary vector perpendicular to v}

followed by normalization. Both formulas need to be used in practice, depending on which one is closer to outputting a degenerate result.