I have two vectors $\vec a$ and $\vec b$ in 3d space. Both of these vectors has length (magnitude) 1 and begin from the origin, so $\vec a$ can be turned into $\vec b$ by a unit quaternion $q$:
$$\vec b = q^{-1} \vec a q$$
So the question is: how can I get this quaternion $q$ that turns vector $\vec a$ into vector $\vec b$? Speaking about why I need this, actually the problem's to rotate a bunch of points, while $\vec a$ and $\vec b$ just set the rotation. And I had to do this exactly using quaternion math, not matrixes or angles.
Any help on how I can solve this would be appreciated, but the better way is to get a rotation quaternion directly without finding a matrix and converting it into a quaternion. Thank in advance!
NOTE: The angle between these two vectors can't be greater than 90°.
I'll use $\hat{a}$ and $\hat{b}$ for the two unit vectors ($\lVert\hat{a}\rVert = 1$ and $\lVert\hat{b}\rVert = 1$. The axis of rotation is then their cross product, and the cosine of the angle is their dot product: $$\begin{aligned} \vec{n} &= \hat{a} \times \hat{b} \\ \hat{n} &= \frac{\vec{n}}{\lVert\vec{n}\rVert} \\ \cos\theta &= \hat{a} \cdot \hat{b} = c, \quad 0 \le \theta \le 180° \\ \cos\frac{\theta}{2} &= \sqrt{\frac{1 + c}{2}} \\ \sin\frac{\theta}{2} &= \sqrt{\frac{1 - c}{2}} \\ \end{aligned}$$ An unit quaternion $\mathbf{q} = (r; i, j, k)$ describes rotation around unit axis vector $\hat{n} = (n_x, n_y, n_z)$ by angle $\theta$, when $$\left\lbrace ~ \begin{aligned} \mathbf{q}_r &= \cos\frac{\theta}{2} = \sqrt{\frac{1 + c}{2}} \\ \mathbf{q}_i &= n_x \sin\frac{\theta}{2} = n_x \sqrt{\frac{1 - c}{2}} \\ \mathbf{q}_j &= n_y \sin\frac{\theta}{2} = n_y \sqrt{\frac{1 - c}{2}} \\ \mathbf{q}_k &= n_z \sin\frac{\theta}{2} = n_z \sqrt{\frac{1 - c}{2}} \\ \end{aligned} \right.$$ As you see, all you need to remember is to normalize the cross product (rotation axis $\vec{n}$) by dividing it by its 2-norm, $\lVert\vec{n}\rVert = \sqrt{\vec{n}\cdot\vec{n}} = \sqrt{n_x^2 + n_y^2 + n_z^2}$), and to halve the rotation angle; I included the two half-angle formulas in the first set for $0 \le \theta \le 180°$.
Also note that the result indeed is an unit quaternion, $\lVert\mathbf{q}\rVert = \sqrt{q_r^2 + q_i^2 + q_j^2 + q_k^2} = 1$. (If weren't isn't, you can always normalize it safely, without introducing any directional bias, by dividing it by its length (that square root). The $q_r = 1$, $q_i = q_j = q_k = 0$ represents no rotation, in case the division yields a non-finite value.)