Representing rotations using quaternions

The key fact is that:

a rotation of an angle $2\theta$ in space, around an axis passing through the origin, is represented by a quaternion $e^{\mathbf{u}\theta}$, where $\mathbf{u}$ is the imaginary quaternion that correspond to the unit vector oriented along the axis of rotation. So we have the correspondence: $$ \vec{w}=R_{\mathbf{u},\theta} \; \vec{v} \quad \longleftrightarrow \quad \mathbf{w}= e^{\mathbf{u}\theta/2}\mathbf{v}e^{-\mathbf{u}\theta/2} $$

See my answer to :Quaternions vs Axis angle.

For the exponentiaton we have that if $ \mathbf{v} \in \mathbb{H}_P$ is an imaginary quaternion, putting $\theta=|\mathbf{v}|$ we have: $$ e^\mathbf{v}= \cos\theta + \mathbf{v}\;\dfrac{\sin \theta}{\theta} $$

See:Exponential Function of Quaternion - Derivation

Why would we need a four dimensional set to represent rotations on a three dimensional space?

Short answer: We don't. A rotation quaternion only has three degrees of freedom, because it is limited to unit length: $$ |Q| = w^2 + x^2 + y^2 + z^2 = 1 $$

Longer answer: If you have any three of those, you can easily calculate the fourth whenever you need it. But if you do that, you lose the possibility of measuring and correcting accumulated errors caused by the limited precision used to represent those numbers in a computer.

After any manipulation of the quaternion, the limited precision may have caused the square sum to deviate somewhat from 1. This deviation gives an estimate of the precision you have lost.

To regain most of that precision, you can multiply each element by the inverse of the length: $$ Q_2 = Q_1 * \sqrt{\frac{1}{w^2 + x^2 + y^2 + z^2}} $$ Thus regaining the unity property by repairing your quaternion as best you can.

Tip for extra high precision: Calculate a re-normalized quaternion after every manipulation. Keep the new one only if its length is closer to 1 than the old one.

Also, don't use float but double.