Assuming that the plane passes through the origin, this problem comes down to finding the intersection of the plane with the nappe of the cone defined by $\mathbf v$ and $\mathbf r$. This intersection will consist of up to two rays. Finding unit vectors in the directions of these rays is easy, from which vectors you can get the corresponding rotation angles or otherwise construct the rotation using your favorite method.
The first thing to do is to check whether or not there is a non-trivial intersection. This can be done by comparing the angles between the rotation axis $\mathbf r$, and $\mathbf v$ and the plane represented by $\mathbf n$, respectively, or, more conveniently, their cosines. Since we’re working with unit vectors, there will be an intersection when $$(\mathbf v\cdot\mathbf r)^2\le1-(\mathbf n\cdot\mathbf r)^2.$$
You can solve the cone-plane intersection problem directly, but I find it easier to transform everything so that the cone is aligned with the $z$-axis. A rotation will do the trick, but I think it’s faster and easier to use a Householder reflection instead. So, compute the transformation matrix $$H=I-2{\mathbf a\mathbf a^T\over\mathbf a^T\mathbf a}$$ where $\mathbf a=\mathbf r-(0,0,1)^T$. This is a reflection in the plane that bisects the angle between $\mathbf v$ and $\mathbf r$. The advantage of using this transformation instead of a rotation is that $H=H^T=H^{-1}=H^{-T}$, so the same matrix can be used to transform back to the original coordinate system and to transform the normal vector $\mathbf n$. (That’s a covariant vector, so if points are transformed as $M\mathbf p$, it transforms as $M^{-T}\mathbf n$.) In practice, you might not need to construct this matrix, since you can directly compute $H\mathbf p=\mathbf p-2{\mathbf p\cdot\mathbf a\over\mathbf a\cdot\mathbf a}\mathbf a$ and only need to map a few vectors.
Set $\mathbf v'=H\mathbf v$ and $\mathbf n'=H\mathbf n$. The vectors that we’re looking for are the intersections of the transformed plane with the circle $x'^2+y'^2=1-(z_{\mathbf v}')^2$, $z'=z_{\mathbf v}'$ (i.e., the circle that is the intersection of the transformed cone nappe with the plane $z'=z_{\mathbf v}'$). Note, by the way, that $z_{\mathbf v}'=\mathbf v\cdot\mathbf r$. Following the method described here, reduce this to a two-dimensional problem by computing the intersection of the transformed target plane with $z'=z_{\mathbf v}'$ and project onto the $x'$-$y'$ plane. For the present problem, this amounts to multiplying the third coordinate of $\mathbf n'$ by $z_{\mathbf v}'$, i.e., $\mathbf l=\operatorname{diag}(1,1,z_{\mathbf v}')\,\mathbf n'$. Proceeding as in the linked answer, find the intersection of this line with the circle $x'^2+y'^2=1-(z_{\mathbf v}')^2$ and project back onto the $z'=z_{\mathbf v}'$ plane to get unit vectors $\mathbf w_1'$ and $\mathbf w_2'$ (which might be identical).
With these vectors in hand, you can proceed in several ways to derive the corresponding rotations. The rotation angle required to take $\mathbf v$ to $\mathbf w_i=H\mathbf w_i'$ is just the negation of the angle between the projections of $\mathbf v'$ and $\mathbf w_i'$ onto the $x'$-$y'$ plane (negated because $H$ is orientation-reversing), and once you have this angle, you can use well-known formulas to construct an appropriate rotation matrix or quaternion. You could also construct a rotation matrix $R'$ in the transformed frame, which will just be a rotation about the $z'$-axis, and then transform it back to the original coordinate system: $R=HR'H$. (This is another instance in which the identity $H=H^T=H^{-1}=H^{-T}$ comes in handy.)
For numerical stability, you might want to align the rotation axis with one of the other coordinate axes. The best choice in that regard is the standard basis vector $\mathbf e_i$ which maximizes $\|\mathbf r-\mathbf e_i\|$. Getting back to the above solution is a simple matter of permuting coordinates after applying $H$, but be sure you keep track of orientation changes if you’re going to use the computed rotation angles directly.
Complex numbers are couples, they have a real part and an imaginay part. Quaternions are hypercomplex numbers, they are also couples of a real part $w$ and an imaginary vector part $b = (x i, y j, z k)$.
Quaternions are tipically parameterized by two parameters: rotation angle $\theta$ and rotation axis (unit vector) $n$. The real part is $w = \cos(\frac{\theta}{2})$ while the vector part is $b = \sin(\frac{\theta}{2}) n$.
Given two qaternions $q_1 = (w_1, b_1)$ and $q_2 = (w_2, b_2)$ their multiplication in vector form is:
$$q_1 q_2 = (w_1, b_1) (w_2,b_2)$$
$$q_1 q_2 = (w_1 w_2 - b_1 \cdot b_2, w_2 b_1 + w_1 b_2 + b_1 \times b_2)$$
After multiplication the new rotation axis is:
$$b_3 = (w_2 b_1 + w_1 b_2 + b_1 \times b_2)$$
$b_3$ has one component in the plane spanned by $b_1$ and $b_2$ ($w_2 b_1 + w_1 b_2$) and one component in the direction orthogonal to that plane ($b_1 \times b_2$).
The case of the SLERP ia a little bit different. Looking at the SLERP formulae:
$$q(t) = \frac{1}{\sin(\frac{\theta}{2})}( \sin((1-t)\frac{\theta}{2}) q_1 + \sin(t \frac{\theta}{2}) q_2)$$
Or
$$q(t) = \frac{1}{\sin(\frac{\theta}{2})}( \sin((1-t)\frac{\theta}{2}) (w_1, b_1) + \sin(t \frac{\theta}{2}) (w_2, b_2))$$
Where $\theta = 2 \cos^{-1}(n_1 \cdot n_2)$ is the angle formed between the axis of rotation of $q_1$ and $q_2$. So the interpolated axis of rotation is:
$$b_t = \frac{1}{\sin(\frac{\theta}{2})}( \sin((1-t)\frac{\theta}{2}) b_1 + \sin(t \frac{\theta}{2}) b_2)$$
Which is in the plane spanned by $b_1$ and $b_2$ (in between $b_1$ and $b_2$ to be precise), the expected behavior of SLERP.
Best Answer
For the equivalence of 1 and 2, we can show that (left or right) translation of the unit quaternions by a unit quaternion is an isometry of $S^3$:
Given vectors $x,y\in S^3$, translating them both by $a\in S^3$ in the sense of quaternion multiplication leads to the dot product
$$ \begin{align} (ax)\cdot(ay) &= \pmatrix{a_0x_0-\vec a\cdot\vec x\\a_0\vec x+x_0\vec a-\vec a\times\vec x} \cdot \pmatrix{a_0y_0-\vec a\cdot\vec y\\a_0\vec y+y_0\vec a-\vec a\times\vec y} \\ &=a_0^2x_0y_0+\vec a^2x_0y_0+a_0^2\vec x\cdot\vec y+(\vec a\cdot\vec x)(\vec a\cdot\vec y)+ (\vec a\times\vec x)(\vec a\times\vec y) \\ &=(a_0^2+\vec a^2)x_0y_0+(a_0^2+\vec a^2)\vec x\cdot\vec y \\ &=x_0y_0+\vec x\cdot\vec y \\ &=x\cdot y\;, \end{align} $$
so the translation is an isometry of $S^3$ and leaves the distribution 2 invariant. (I'm taking the more technical aspects of the Haar measure involving inner and outer regularity for granted in this context.)
If we split the rotations up as in 3, then for the distribution to be invariant under translation by arbitrary rotations, in particular the distribution of $\theta$ must be invariant under arbitrary rotations about $\hat w$, so $\theta$ must be uniformly distributed. Also the distribution of $\hat w$ clearly has to be rotationally symmetric about $\hat z$. Thus, it remains only to show that the uniform distribution of $\hat w\cdot\hat z$ in distribution 3 is correct.
I don't have an elegant argument without calculus for this, but here's a calculation that tries to do without too much brute force. For a rotation by $\alpha$ about $\hat u$, the rotated $z$ axis is
$$ \cos\alpha\,\hat z+\sin\alpha\,\hat u\times\hat z+(1-\cos\alpha)\hat u(\hat u\cdot\hat z)\;, $$
and forming the dot product with $\hat z$ yields a relationship between the three cosines,
$$c=a+(1-a)b^2\;,$$
where $a=\cos\alpha$, $b=\hat u\cdot\hat z$ and $c$ is the cosine between the original and rotated $z$ axes.
To calculate the cumulative distribution function for $c$, note that $c\ge a$ and that $b$ is uniformly distributed according to distribution 2, so we can solve for the upper limit of $b$,
$$ b=\sqrt{\frac{c-a}{1-a}}\;, $$
and find the density for $a$ according to distribution 2,
$$ \sin^2\frac\alpha2\mathrm d\alpha\propto\frac{1-\cos\alpha}{\sin\alpha}\mathrm d(\cos\alpha)=\frac{1-a}{\sqrt{1-a^2}}\mathrm da=\sqrt{\frac{1-a}{1+a}}\mathrm da\;, $$
to find the cumulative distribution function for $c$ proportional to
$$ \int_{-1}^c\sqrt{\frac{c-a}{1-a}}\sqrt{\frac{1-a}{1+a}}\mathrm da=\int_{-1}^c\sqrt{\frac{c-a}{1+a}}\mathrm da\;. $$
Now we can differentiate with respect to $c$ to get the density for $c$. The term from varying the integration limit vanishes since the integrand vanishes at $a=c$, so the result is proportional to
$$ \int_{-1}^c\frac1{\sqrt{(c-a)(1+a)}}\mathrm da=\left[2\arctan\sqrt{\frac{1+a}{c-a}}\right]_{-1}^c=\pi\;, $$
a constant, which establishes the equivalence between distributions 2 and 3.