The simplest way to model 3D rotations with quaternions is to view 3D space as the "imaginary quaternions", that is, use $i,j,k$ vectors like you do in physics to model 3D space. (But no real part.)
Every quaternion with norm 1 produces a rotation of this 3D space by conjugation $x\mapsto qxq^{-1}$, where $q$ is any length 1 quaternion.
The formula is surprisingly simple: pick a vector that points along the axis you want to rotate around, and scale it down to a unit vector, call it $h$. Then, take your angle of rotation $\theta$ and form $q=\cos(\theta/2)+h\sin(\theta/2)$. The result is a unit length quaternion that produces the rotation you wanted. ($\theta$ may be positive or negative, depending on the direction of rotation.)
I'm not sure I understand what you mean by "nullify rotation in the x and y axes", but I was guessing you meant you would like to be able to pick which axis you are rotating around.
This is called the decomposition into Euler angles. The problem is that the result depends on which order you're planning to do the x-, y-, and z-rotations in.
Here's the drill, though:
Step 1: Convert your rotation to a $3 \times 3$ orthogonal matrix using Rodrigues' formula. This is well-detailed on Wikipedia and elsewhere, so I won't go into it here.
Step 2: Your matrix $M$ has three columns, $u_1, u_2, u_3$. Let's say you want to write
$$
M = R_x R_y R_z
$$
where the $R$s are rotations (by varying amounts) around the three axes. We could then write
$$
R_y^{-1}R_x^{-1}M = R_z
$$
which, since a rotation about $z$ has interesting entries only in its upper left $2 \times 2$ block, means that we need
$$
R_y^{-1}R_x^{-1}u_3 = e_3
$$
where $e_3$ is the vector $\begin{bmatrix}0\\0\\1\end{bmatrix}$.
Let's make that happen. First, we'll find an $x$-rotation that "kills off" the first entry of $u_3$.
Case 1: Both the first and last entries of $u_3$ are 0. In this case, perform no rotation about the $y$ axis, and if the second entry is 1, pick $R_x$ to rotate about $x$ by $90$ degrees; if the second entry is $-1$, pick $R_x$ to rotate about $x$ by $-90$ degrees. With these choices, $R_x^{-1}u_3$ will be $e_3$.
Case 2: At least one of these entries is nonzero. Let's call them $s$ and $c$. Let $\alpha = \text{atan2}(c, s)$, and let $R_y$ be rotation about $y$ by angle $\alpha$. Then $R_y^{-1} u_3$ will have its first entry zero, and its third entry nonzero. Call the second and third entries of $R_y^{-1} u_3$ by the names $a$ and $b$ respectively. Let $\beta = \text{atan2}(b, a)$, and let $R_x$ be rotation about $x$ by angle $\beta$. Then
$$
R_x^{-1}R_y^{-1} u_3 = e_3.
$$
That means that
$$
R_x^{-1}R_y^{-1} M = \begin{bmatrix}
c' & -s' & 0\\
s' & c' & 0\\
0 & 0 & 1
\end{bmatrix},
$$
and we can write $\gamma = \text{atan2}(c', s')$ and let $R_z$ be rotation about $z$ by angle $\gamma$, and we're done.
It's quite possible that I've made a sign error here somewhere: that one of the angles should be the negative of what I wrote. But if you write a program to implement this, it should be pretty obvious where any error is.
Best Answer
The angle between two unit vectors $u$ and $v$ is given as
$$ \theta = \arccos\left(u^\top v\right) $$
the shortest axis of rotation is the vector orthogonal to both vectors, and can be found using the cross product
$$ s = \frac{u \times v}{\lVert u \times v\rVert} $$
In most cases, $s$ should be normalized, even if $u$ and $v$ are already both unit vectors.
(note that switching the order of the cross product will give you the negative of the axis. The right-hand rule determines the order)
To create the rotation matrix that rotates $u$ into $v$, we can use the axis-angle formula for rotation matrices: (taken from wikipedia)
$$ R = \begin{bmatrix} \cos\theta + s_x^2 (1-\cos\theta) & s_x s_y (1-\cos\theta) -s_z \sin\theta & s_x s_z(1-\cos\theta + s_y \sin \theta \\ s_y s_x (1-\cos\theta)+s_z \sin\theta & \cos\theta + s_y^2 (1-\cos\theta) & s_y s_z (1-\cos\theta) - s_x\sin\theta \\ s_z s_x (1-\cos\theta) - s_y \sin\theta & s_zs_y(1-\cos\theta) + s_x \sin\theta & \cos\theta + s_z^2(1-\cos\theta) \end{bmatrix} $$
and now you can compute
$$ v = Ru $$
and
$$ u = R^\top v $$