>>16773205
Possible solution?
Start from the axis of revolution [math] \mathbf{a} = 0 + a_x \mathbf{i} + a_y \mathbf{j} + a_z \mathbf{k} [/math], [math] \| a \| = 1 [/math]
and the CCW angle [math] \theta [/math] bounded between [math] 0 [/math] and [math] 2\pi [/math].
The Wikipedia formula works and is numerically stable, but for small thetas you can't extract the axis from the quaternion. An angle of 0.25 degrees or less will give you NaN for the axis, based on my Python experiment.
[eqn]
\begin{align}
\mathbf{p} &= p_r + p_i \mathbf{i} + p_j \mathbf{j} + p_k \mathbf{k} \\
p_r &= \cos(\frac{1}{2}\theta) \\
p_i &= a_x \sqrt{1 - p_r^2} \\
p_j &= a_y \sqrt{1 - p_r^2} \\
p_j &= a_z \sqrt{1 - p_r^2}
\end{align}
[/eqn]
Alternatively: let's just define [math] m \in \mathbf{R} [/math] such that [math] \| \mathbf{p} - p_r \| = m \|\mathbf{a}\| [/math]
ie. [math] m = \frac{ \| \mathbf{p} - p_r \| }{ \|\mathbf{a}\| } = \frac{\sqrt{p_i^2 + p_j^2 + p_k^2}}{1} = \sqrt{1 - p_r^2} [/math].
Rewrite [math] \mathbf{p} [/math] like this:
[eqn]
\begin{align}
\mathbf{p} &= \frac{m}{m} p_r + m a_x \mathbf{i} + m a_y \mathbf{j} + m a_z \mathbf{k} \\
&= m \left( \frac{p_r}{m} + a_x \mathbf{i} + a_y \mathbf{j} + a_z \mathbf{k} \right)
\end{align}
[/eqn]
Set
[eqn]
\mathbf{q} = \frac{p_r}{m} + a_x \mathbf{i} + a_y \mathbf{j} + a_z \mathbf{k}
[/eqn]
Rotation formula:
[eqn]
\begin{align}
R_\mathbf{p}(\mathbf{v}) &= \mathbf{p} \mathbf{v} \mathbf{p}^* \\
&= (m \mathbf{q})\ \mathbf{v}\ (m \mathbf{q})^* \\
&= (m \mathbf{q})\ \mathbf{v}\ (\mathbf{q}^* m^*) \\
&= m\ (\mathbf{q} \mathbf{v} \mathbf{q}^*)\ m \\
&= (m^2)(\mathbf{q} \mathbf{v} \mathbf{q}^*)
\end{align}
[/eqn]
Where [math] m^2 = 1 - p_r^2 [/math].
If I do this I will always have access to the axis of revolution in exchange for the more expensive real term [math] \frac{p_r}{m} = \frac{p_r}{\sqrt{1 - p_r^2}} [/math].
