The correct way to handle quaternions and rotations? Let see if I should post it here or go to Stack Overflow.
I'm trying to rotate a point $p_0$ on a unit sphere that depends on time. Which is the "authentic/correct" way to do so? Use the spatial rotation algorithm: $$\text{Rotate}(p_0, q_1, t) = (q_1)^t \otimes p_0 \otimes (q_1)^{-t} \text{ } (\text{Rotate}(p_0, q_1, 1) = q_1 \otimes p_0 \otimes (q_1)^{-1} = p_1),$$ or the SLERP algorithm: $$\text{SLERP}(p_1, p_0, t) = (p_0) \otimes ((p_0)^{-1} p_1)^t?$$
 A: Your $\sf Rotate$ algorithm seems to rotate by $2t$, not $t$. It ought to be
$$ \mathsf{Rotate}(t)=q^{t/2}p_0q^{-t/2}. $$
Here, the initial and final positions are $p_0$ and $p_1$ respectively, $q=p_1p_0^{-1}=e^{\Omega\mathbf{v}}=\cos(\Omega)+\sin(\Omega)\mathbf{v}$ with $\mathbf{v}$ a rotation axis (i.e. a 3D unit vector), and $q^{t/2}:=e^{(t/2)\Omega\mathbf{v}}$. As $t$ varies in the interval $[0,1]$, this will trace out the arc from $p_0$ to $p_1$ on the 2D sphere $S^2$ within 3D Euclidean space $\mathbb{R}^3$. Note it matches 3D SLERP:
$$ \mathsf{Rotate}(t)=\frac{\sin((1-t)\Omega)}{\sin \Omega}p_0+\frac{\sin(t\Omega)}{\sin\Omega}p_1. $$
The formulas must match because both describe the unique point a fraction of $t$ along the arc from $p_0$ to $p_1$ (assuming, of course, $p_1\ne \pm p_0$ is distinct and not antipodal).
Note that, because $p_0$ and $p_1$ are unit vectors and the imaginary part of $q=p_1p_0^{-1}$ is their cross product, the imaginary part $\mathbf{v}=p_0\times p_1$ is perpendicular to $p_0$ and $p_1$, which means $\mathbf{v}$ anticommutes with them, so we can slide $q$ (or powers of $q$) across $p_0$ or $p_1$ if we invert $q$. So we have $q^tp_0=q^{t/2}p_0q^{-t/2}=p_0q^{-t}$.
Therefore, the SLERP algorithm $\mathsf{slerp}(t)=p_0q^{-t}=p_0(p_0^{-1}p_1)^t$ (note because $p_0$ and $p_1$ are vectors, $q^{-1}=p_0p_1^{-1}=-p_0p_1=p_0^{-1}p_1$) trace the same arc in $S^2$. But evidently $q^tp_0$ and $p_0q^{-t}$ require less multiplication than the conjugation $q^{t/2}p_0q^{-t/2}$ so it would seem (to my untrained eye) to be less computationally expensive to use the SLERP algorithm than the Rotate algorithm.
Another thing to keep in mind is that, for given $q$, the function $f(p)=q^{t/2}pq^{-t/2}$ rotates all unit vectors $p$ by the angle $t\Omega$ around $\mathbf{v}$, however the same is not true of $g(p)=q^tp$ or $h(p)=pq^t$; the latter functions trace out trajectories in the three-sphere $S^3$ in 4D Euclidean space (the quaternions $\mathbb{H}$), not $S^2$. That is, for general $p$ and $q$, the expressions $q^tp$ and $pq^t$ are not 3D vectors, although they are unit quaternions.
