An unit quaternion $\mathbf{q}$ that rotates $v$ to $v^\prime$ (given $\lVert v \rVert = \lVert v^\prime \rVert$) can be defined as
$$\mathbf{q} = \cos\frac{\theta}{2} + a \sin\frac{\theta}{2}$$
where
$$\begin{aligned}
a &= \frac{v \times v^\prime}{\left\lVert v \times v^\prime \right\rVert} \\
\cos \theta &= \frac{v \cdot v^\prime}{\left\lVert v \right\rVert \left\lVert v^\prime \right\rVert} \\
\sin \theta &= \frac{\left\lVert v \times v^\prime \right\rVert}{\left\lVert v \right\rVert \left\lVert v^\prime \right\rVert} \\
\end{aligned}$$
To additionally rotate $u$ to $u^\prime$, you need a rotation $\mathbf{p}$ around axis $v^\prime$ by suitable angle $\varphi$.
Let $u_q = \mathbf{q} u \mathbf{q}^{-1}$ (where $\mathbf{q}^{-1} = \sin\frac{\theta}{2} - a \cos\frac{\theta}{2}$), i.e. $u$ rotated by $\mathbf{q}$. Then,
$$\begin{aligned}
b &= \frac{v^\prime}{\left\lVert v^\prime \right\rVert} \\
\cos\varphi &= \frac{u_q \cdot u^\prime}{\left\lVert u_q \right\rVert \left\lVert u^\prime \right\rVert} \\
\mathbf{p} &= \cos\left(\frac{\varphi}{2}\right) + b \sin\left(\frac{\varphi}{2}\right) \\
\end{aligned}$$
and the combined rotation needed is $\mathbf{p}\mathbf{q}$.
Another, more general option is to construct the two orientations' basis vectors (rotation matrices), combine them to get the rotation matrix needed, then recover the rotation quaternion from the rotation matrix.
If you have two linearly independent vectors $\vec{a}$ and $\vec{b}$ (meaning there is no $\lambda \in \mathbb{R}$ so that $\vec{a} = \lambda \vec{b}$), you can construct the basis vectors $\hat{e}_1$, $\hat{e}_2$, and $\hat{e}_3$ very easily. The first basis vector is just one of the vectors scaled to unit length:
$$\hat{e}_1 = \frac{\vec{a}}{\left\lVert\vec{a}\right\rVert}$$
The second basis vector is the perpendicular part of the second vector with respect to the first basis vector. We can obtain this by one step of the Gram–Schmidt process, and normalizing the result to unit length:
$$\hat{e}_2 = \frac{ \vec{b} - \vec{e}_1 \left( \vec{e}_1 \cdot \vec{b} \right) }{ \left\lVert \vec{b} - \vec{e}_1 \left( \vec{e}_1 \cdot \vec{b} \right) \right\rVert }$$
The third basis vector is the cross product of the two:
$$\hat{e}_3 = \hat{e}_1 \times \hat{e}_2$$
and the rotation matrix describing that is
$$\mathbf{R} = \left[ \begin{matrix} \hat{e}_1 & \hat{e}_2 & \hat{e}_3 \end{matrix} \right] = \left[ \begin{matrix}
x_1 & x_2 & x_3 \\
y_1 & y_2 & y_3 \\
z_1 & z_2 & z_3 \\
\end{matrix} \right ]$$
Because this matrix is orthonormal, its inverse is its transpose, $\mathbf{R}^{-1} = \mathbf{R}^T$.
If $\mathbf{R}_1$ describes the current orientation, and $\mathbf{R}_2$ the desired orientation, then
$$\mathbf{R}_{1\to 2} = \mathbf{R}_2 \mathbf{R}_1^{-1} = \mathbf{R}_2 \mathbf{R}_1^T$$
is the rotation needed.
To recover quaternion $$\mathbf{q} = w + x\mathbf{i} + y\mathbf{j} + z\mathbf{k}$$ from a pure 3×3 rotation matrix
$$\mathbf{R} = \left[ \begin{matrix}
r_{11} & r_{12} & r_{13} \\
r_{21} & r_{22} & r_{23} \\
r_{31} & r_{32} & r_{33} \\
\end{matrix} \right]$$
in a robust, numerically stable manner, I recommend the following procedure:
If $r_{11} + r_{22} + r_{33} \ge 0$:
$$\left\lbrace ~ \begin{aligned}
w &= \sqrt{1 + r_{11} + r_{22} + r_{33}} / 2\\
x &= \pm \sqrt{1 + r_{11} - r_{22} - r_{33}} / 2, \text{ same sign as } r_{32} - r_{23} \\
y &= \pm \sqrt{1 - r_{11} + r_{22} - r_{33}} / 2, \text{ same sign as } r_{13} - r_{31} \\
z &= \pm \sqrt{1 - r_{11} - r_{22} + r_{33}} / 2, \text{ same sign as } r_{21} - r_{12} \\
\end{aligned} \right.$$
Otherwise, if $r_{11} \ge r_{22}$ and $r_{11} \ge r_{33}$:
$$\left\lbrace ~ \begin{aligned}
s &= 2 \sqrt{ 1 + r_{11} - r_{22} - r_{33}} \\
w &= ( r_{32} - r_{23} ) / s \\
x &= s / 4 \\
y &= ( r_{21} + r_{12} ) / s \\
z &= ( r_{13} + r_{31} ) / s \\
\end{aligned} \right.$$
Otherwise, if $r_{22} \ge r_{11}$ and $r_{22} \ge r_{33}$:
$$\left\lbrace ~ \begin{aligned}
s &= 2 \sqrt{ 1 - r_{11} + r_{22} - r_{33}} \\
w &= ( r_{13} - r_{31} ) / s \\
x &= ( r_{21} + r_{12} ) / s \\
y &= s / 4 \\
z &= ( r_{32} + r_{23} ) / s \\
\end{aligned} \right.$$
Otherwise:
$$\left\lbrace ~ \begin{aligned}
s &= 2 \sqrt{ 1 - r_{11} - r_{22} + r_{33}} \\
w &= ( r_{21} - r_{12} ) / s \\
x &= ( r_{13} + r_{31} ) / s \\
y &= ( r_{32} + r_{23} ) / s \\
z &= s / 4 \\
\end{aligned} \right.$$
This is based on the quaternion-derived rotation matrix $\mathbf{R}$,
$$\mathbf{R} = \left[ \begin{matrix}
1 - 2 (y^2 + z^2) & 2 (x y - w z) & 2 (x z + w y) \\
2 (x y + w z) & 1 - 2 (x^2 + z^2) & 2 (y z - w x) \\
2 (x z - w y) & 2 (y z + w x) & 1 - 2 (x^2 + y^2) \\
\end{matrix} \right] = \left[ \begin{matrix}
r_{11} & r_{12} & r_{13} \\
r_{21} & r_{22} & r_{23} \\
r_{31} & r_{32} & r_{33} \\
\end{matrix} \right]$$
having the following properties:
$$\begin{aligned}
r_{11} &= w^2 + x^2 - y^2 - z^2 \\
r_{22} &= w^2 - x^2 + y^2 - z^2 \\
r_{33} &= w^2 - x^2 - y^2 + z^2 \\
r_{32} - r_{23} &= 4 w x \\
r_{13} - r_{31} &= 4 w y \\
r_{21} - r_{12} &= 4 w z \\
r_{21} + r_{12} &= 4 x y \\
r_{13} + r_{31} &= 4 x z \\
r_{32} + r_{23} &= 4 y z \\
r_{11} + r_{22} + r_{33} &= 3 w^2 - x^2 - y^2 - z^2 \\
1 + r_{11} + r_{22} + r_{33} &= 4 w^2 ~ \text{ if } ~ w = \sqrt{1 - x^2 - y^2 - z^2} \\
1 + r_{11} - r_{22} - r_{33} &= 4 x^2 ~ \text{ if } ~ w = \sqrt{1 - x^2 - y^2 - z^2} \\
1 - r_{11} + r_{22} - r_{33} &= 4 y^2 ~ \text{ if } ~ w = \sqrt{1 - x^2 - y^2 - z^2} \\
1 - r_{11} - r_{22} + r_{33} &= 4 y^2 ~ \text{ if } ~ w = \sqrt{1 - x^2 - y^2 - z^2} \\
\end{aligned}$$
When $r_{11} + r_{12} + r_{13} \ge 0$, we can recover the quaternion from the diagonal entries relying on $w^2 + x^2 + y^2 + z^2 = 1$. However, when the sum is negative $w$ is very small, and we can get better numerical stability by starting with $x^2$, $y^2$, or $z^2$ (whichever is greatest), and using the non-diagonal elements of $\mathbf{R}$ to obtain the other three components.
