How to decompose a unit quaternion into 3 Tait-Bryan quaternions instead of 3 real numbers?

I'm familiar with the formulas for decomposing a unit quaternion $Q$ into chained Tait-Bryan angles $\phi\theta\psi$ (Wikipedia has the formulas for the $zyx$ chain here), but I'm looking to instead decompose $Q$ into three quaternions $A \times B \times C$, each representing a rotation around a single axis (and therefore composed of only two nonzero terms, the real coefficient and one of three imaginary coefficients). $ABC$ can be reverse-engineered from Tait-Bryan angles from the formulas above, but this has two problems:

1. (EDIT: This is incorrect, but I'm not sure how.) These formulas can't represent $Q_w < -\frac{\sqrt2}{2}$ without inverting $Q$ (i.e. $A \times B \times C = -Q$). The formulas output $|\phi| \le \pi$, $|\theta| \le \frac{\pi}{2}$, $|\psi| \le \pi$, but when $Q_w < -\frac{\sqrt2}{2}$, the rotation chain can't be represented within those bounds without inversion.
2. Reverse-engineering $ABC$ requires what seems to be superfluous trig (taking the sines & cosines of halved arcsines or arctangents). This should only require matrix operations and not any trig, right?

I believe that there exist 8 unique $ABC$ for any given $Q$ (they involve inversion and reversing terms of $ABC$), but I haven't worked on a concise, digestable representation of these. However, for all $Q$ for $Q_w > -\frac{\sqrt2}{2}$, exactly one of these $ABC$ represents $\phi\theta\psi$ within the bounds described above.

I do also know how to represent $Q$ in terms of $ABC$:

Chained rotation $A \times B \times C = Q$ about axes $tuv$, for $tuv = xyz, yzx,$ or $zxy$:

$A_wB_wC_w - A_tB_uC_v = Q_w$

$A_tB_wC_w + A_wB_uC_v = Q_t$

$A_wB_uC_w - A_tB_wC_v = Q_u$

$A_wB_wC_v + A_tB_uC_w = Q_v$

(EDIT: Expression with B factored out: $\frac{A_wC_w - A_tC_v}{A_tC_w + A_wC_v} = \frac{Q_w + Q_u}{Q_t + Q_v}$)

Chained rotation $A \times B \times C = Q$ about axes $tuv$, for $tuv = xzy, yxz,$ or $zyx$ (same as above with flipped signs):

$A_wB_wC_w + A_tB_uC_v = Q_w$

$A_tB_wC_w - A_wB_uC_v = Q_t$

$A_wB_uC_w + A_tB_wC_v = Q_u$

$A_wB_wC_v - A_tB_uC_w = Q_v$

But I'm stuck here. Can someone point me in the right direction?

EDIT: I'm also aware that it's still impossible to differentiate between $A$ and $C$ at the Euler singularities (i.e. due to gimbal lock).

Okay, this is what I wound up with. This replicates the results of the Tait-Bryan decomposition formulas without any trig. Chained rotation $A×B×C=Q$ about axes $tuv,$ for $tuv=xyz,yzx,$ or $zxy$:

$B_u = sgn(Q_wQ_u + Q_tQ_v) \times \sqrt{\frac{1 - \sqrt{1 - 4(Q_wQ_u + Q_tQ_v)^2}}{2}}$

$B_w = \sqrt{1 - B_u^2}$

$A_t = sgn(Q_wQ_t - Q_vQ_u) \times \sqrt{\frac{\frac12 - Q_u^2 - Q_t^2}{1 - 2B_u^2} + \frac12}$

$A_w = \sqrt{1 - A_t^2}$

$C_v = sgn(Q_wQ_v - Q_tQ_u) \times \sqrt{\frac{\frac12 - Q_u^2 - Q_v^2}{1 - 2B_u^2} + \frac12}$

$C_w = \sqrt{1 - C_v^2}$