Let $\mathbf{a}$, $\mathbf{b}$, $\mathbf{c}$ be points on the unitary sphere, so that $\|\mathbf{a}\| = \| \mathbf{b} \| = \| \mathbf{c} \| = 1$.
Let $\mathbf{a'}$, $\mathbf{b'}$, $\mathbf{c'}$ be points on the same sphere, so that:
$\mathbf{a}\cdot\mathbf{b} = \mathbf{a'} \cdot \mathbf{b'}$
$\mathbf{a}\cdot\mathbf{c} = \mathbf{a'} \cdot \mathbf{c'}$
$\mathbf{b}\cdot\mathbf{c} = \mathbf{b'} \cdot \mathbf{c'}$
Therefore, a rotation exists so that the first three versors can be mapped on the other three. In general a rotation in $\mathbb{R}^3$ can be expressed as
$R_{\alpha\beta\gamma} = R_x(\gamma) R_y(\beta) R_z(\alpha)$
Where $\alpha$, $\beta$ and $\gamma$ are the Euler angles for the rotation (yaw, pitch and roll). Now, I need to determine $\alpha$, $\beta$ and $\gamma$ so that:
$\mathbf{a} = R_{\alpha\beta\gamma} \mathbf{a'}$
$\mathbf{b} = R_{\alpha\beta\gamma} \mathbf{b'}$
$\mathbf{c} = R_{\alpha\beta\gamma} \mathbf{c'}$
I tried to carry out the explicit computation but it turns out that it is quite a pain to do it. Is there a smart way to do such a thing?
Setting aside issues of representation for the moment, the most straightforward way to solve this is by finding the rotation $\mathcal{P}$ that takes the initial $\mathbf{a}$ and $\mathbf{b}$ to 'canonical' locations, and similarly finding the rotation $\mathcal{Q}$ that takes $\mathbf{a'}$ and $\mathbf{b'}$ to those same canonical axes; the rotation that takes $\mathbf{a}$ to $\mathbf{a'}$ and $\mathbf{b}$ to $\mathbf{b'}$ can then be written as $\mathcal{R} = \mathcal{P}\cdot \mathcal{Q^{-1}}$ (or possibly as $\mathcal{Q^{-1}}\cdot\mathcal{P}$, depending on how you choose to take rotations 'acting' — this is a slightly tricky issue with no canonical answer!).
Finding $\mathcal{P}$ (and by extension $\mathcal{Q}$) similarly breaks down into two distinct steps: (a) find a rotation $\mathcal{P}_{\mathbf{a}}$ that takes $\mathbf{a}$ to (e.g.) the $x$-axis, and then (b) find a rotation $\mathcal{P}_{\mathbf{b}}$ that takes $\mathbf{b}$ to lie on the $xy$ plane (and canonically, in the $+y$ half-plane).
The first of these can be done by finding the vector $\mathbf{v_a} = \mathbf{a}\times \mathbf{e}_x$ (where $\mathbf{e}_x = (1,0,0)$ is the unit $x$ vector); the direction of this vector gives an axis to rotate around, while its magnitude helps determine the angle of rotation (according to the classic $\left|\mathbf{u}\times \mathbf{v}\right| = \left|\mathbf{u}\right|\left|\mathbf{v}\right|\sin(\theta)$ formula). There's one small ambiguity, though: this formula doesn't allow for distinguishing $\theta$ from $\pi-\theta$, since the sine of both values is the same. To clear that up, look at the sign of $\mathbf{a}\cdot \mathbf{e}_x$ (i.e., of $a_x$) to determine whether the angle between $\mathbf{a}$ and the $x$ axis is acute or obtuse.
To find the second, first apply the rotation that was just found to the vector $\mathbf{b}$ to generate a new vector $\mathbf{d} = \mathcal{P}_{\mathbf{a}}(\mathbf{b})$; we then need to find the amount of rotation about the $x$-axis that brings $\mathbf{d}$ into the positive $xy$-plane. But of course, this rotation is just $\phi = \mathrm{atan2}(d_z, d_y)$, and so we get that $\mathcal{P}_{\mathbf{b}}$ is the rotation about the $x$ axis by this amount. A similar approach finds the rotation $\mathcal{Q}$ from $\mathbb{a'}$ and $\mathbb{b'}$.
Note that when you implement this, you don't want to actually get the angles via asin/atan; instead, use trig identities to express the sine and cosine values that you need to calculate for the quaternions in terms of the quantities you're taking the arcsin of (and then adjust for signs, as I mentioned above).