Decomposing a quarternion into unit axes rotations

As has already been written in the comments, there are singular points at which this inverse is ill-defined, and the point of using quaternions is to avoid these sorts of problems. In case you can't avoid this inversion (e.g. because you're using other people's code), here's an idea for how you might go about it:

With $\alpha=X/2$, $\beta=Y/2$, $\gamma=Z/2$, the quaternion $Q$ comes out as

$$ \begin{pmatrix} \cos\alpha\cos\beta\cos\gamma-\sin\alpha\sin\beta\sin\gamma\\ \sin\alpha\cos\beta\cos\gamma+\cos\alpha\sin\beta\sin\gamma\\ \cos\alpha\sin\beta\cos\gamma-\sin\alpha\cos\beta\sin\gamma\\ \sin\alpha\sin\beta\cos\gamma+\cos\alpha\cos\beta\sin\gamma\\ \end{pmatrix} $$

Forming squares and combining them in various ways to simplify the trigonometric functions yields

$$ \begin{eqnarray} q_0^2-q_1^2-q_2^2+q_3^2&=&\cos X\cos Y,\\ q_0^2+q_1^2-q_2^2-q_3^2&=&\cos Y\cos Z,\\ q_0^2-q_1^2+q_2^2-q_3^2&=&\cos X\cos Z-\sin X\sin Y\sin Z . \end{eqnarray} $$

If you take $\sin X\sin Y\sin Z$ to one side in the last equation and square, you can express the squared sines by squared cosines and then substitute $\cos X$ and $\cos Z$ from the first two equations. This yields a cubic equation for $\cos^2 Y$. You'll still have to deal with all the usual issues about the signs and ranges of the angles, but this should at least solve the algebraic part of the problem. (Regarding the problem of equivalent angles, see also How can I find equivalent Euler angles?)