In $\mathbb R^4$ a reflect in a hyperplane through the origin $O$ is, $\forall q \in \mathbb R^4$ $$q \mapsto -u\overline{q}u$$ , where $u$ is a unit quaternion.
In $\mathbb R^3$ a rotation fixing the origin $O$ is $$p \mapsto x^{-1}px$$ , where $p$ is a pure quaternion, and $x$ is a unit quaternion.
Since a rotation is equivalent to 2 reflects chained up, given unit quaternion $u$ for a rotation $p \mapsto u^{-1}pu$, I'm trying to find unit quaternions $x$, $y$ so that the rotation can be decomposed to two reflects, first via $x$, then via $y$, i.e.
$$p \mapsto -x\overline{p}x \mapsto -y \overline{-x\overline{p}x} y = yx^{-1}px^{-1}y$$
Let this equals to $u^{-1}pu$ we get
$$x^{-1}y=u, \qquad yx^{-1}=u^{-1}$$
Please enlighten me how to find such a pair of $x$ and $y$ ?
In 2D the rotation by an angle of $\theta$ is achieved by reflecting across two lines (intersecting at the center of rotation) separated by an angle of $\frac{1}{2}\theta$. (The order in which you do the reflections determines the orientation of the rotation.) This generalizes to 3D with two planes intersecting at the axis of rotation with a dihedral angle of $\frac{1}{2}\theta$, and it generalizes to any number of dimensions with hyperplanes (i.e. codimension $1$ subspaces) intersecting at the orthogonal complement of a 2D plane of rotation.
Every unit quaternion, aka versor, has a polar form $e^{\theta\mathbf{u}}$ with $\mathbf{u}\in S^2\subset\mathbb{R}^3$ and $0\le\theta\le\pi$ is convex, which is unique except for real versors $\pm1$ (corresponding to $\theta=0,\pi$ respectively and with $\mathbf{u}\in S^2$ arbitrary), and by Euler equals $\cos(\theta)+\sin(\theta)\mathbf{u}$. (Note the square roots of $-1$ in the quaternions are precisely the unit vectors.) Conjugating a 3D vector by $e^{\theta\mathbf{u}}$ has the effect of rotating it around the $\mathbf{u}$-axis by $2\theta$. To see this, describe everything wrt an orthonormal basis $\{1,\mathbf{u},\mathbf{v},\mathbf{w}\}$.
When $p$ is a versor we can show $x\mapsto -p\overline{x}p$ is a reflection across the 3D hyperplane perpendicular to $p$. (We can get this formula by seeing $-\overline{x}$ works for $p=1$ and then conjugating this map by either left- or right-multiplication by $p$.) So to rotate 3D space around an axis $\mathbf{u}$ by an angle $\theta$, we can pick two vectors $\mathbf{v}$ and $\mathbf{w}$ which are orthogonal to $\mathbf{u}$ but at an (correctly oriented) angle of $\frac{1}{2}\theta$ to each other, then reflect across $\mathbf{u}^\perp$ followed by reflecting across $\mathbf{v}^\perp$.