Suppose we want to calculate the cross product $u\times v$ for two vectors $u,v\in\mathbb{R}^3$ on a computer satisfying $x\circledast y=(x*y)(1+\epsilon)$, with $|\epsilon|\leq\epsilon_{\mathrm{machine}}$ for $*$ one of the operations $\cdot,+,-$; and $\mathrm{fl}(x)=x(1+\epsilon)$ with $|\epsilon|\leq\epsilon_{\mathrm{machine}}$.
A backward stable algorithm for this problem can be found as follows. Write $u=\begin{bmatrix}u_1&u_2&u_3\end{bmatrix}^T$. Then define $$f_1:[u|v]\mapsto\begin{bmatrix}\phantom{-}0\\ -\mathrm{fl}(u_1)\otimes\mathrm{fl}(v_3)\\ \phantom{-}\mathrm{fl}(u_1)\otimes\mathrm{fl}(v_2)\end{bmatrix},$$ thus applying the most straightforward cross product algorithm to $\begin{bmatrix}u_1&0&0\end{bmatrix}^T$ and $v$; define $f_2,f_3$ similarly. Then using the linearity of the cross product and assigning the errors to different elements of the $3\times2$ matrix $[u|v]$ for any of the three algorithms, we obtain the backward stable algorithm $\tilde{f}=(f_1\oplus f_2)\oplus f_3$.
The idea is to try to assign the errors of $f_i$ to $u_i$ and $v_{i+1}$ only (indices modulo $3$); we look for some for some $W=[w_1|w_2]\in\mathbb{R}^{3\times2}$, such that $$ \tilde{f}([u|v])=f(W):=w_1\times w_2 ~~\text{ and }~~ \frac{\|[u|v]-W\|}{\|[u|v]\|}=\mathcal{O}(\epsilon_{\mathrm{machine}}). $$ For example, in the second element of $f_1$, we can find $\epsilon_i$ with $|\epsilon_i|\leq\epsilon_{\mathrm{machine}}, i=1,2,3,$ such that $-\mathrm{fl}(u_1)\otimes\mathrm{fl}(v_3)=u_1v_3(1+\epsilon_1)(1+\epsilon_2)(1+\epsilon_3)$; assign these errors to $u_1$ by taking $$w_{11}=u_1(1+\epsilon_1)(1+\epsilon_2)(1+\epsilon_3)=u_1(1+\mathcal{O}(\epsilon_{\mathrm{machine}}))$$ (this will one more time adjusted later). Similarly, in the third element of $f_3$, we can find $\epsilon_i$ with $|\epsilon_i|\leq\epsilon_{\mathrm{machine}}, i=4,5,6,$ such that $$\mathrm{fl}(u_1)\otimes\mathrm{fl}(v_2)=u_1v_2(1+\epsilon_4)(1+\epsilon_5)(1+\epsilon_6);$$ Now we want to choose $w_{22}=\frac{(1+\epsilon_4)(1+\epsilon_5)(1+\epsilon_6)}{(1+\epsilon_1)(1+\epsilon_2)(1+\epsilon_3)}=v_2(1+\mathcal{O}(\epsilon_{\mathrm{machine}}))$ and then $w_{11}w_{22}=\mathrm{fl}(u_1)\otimes\mathrm{fl}(v_2)$ exactly.
We do this similarly for the other algorithms $f_1$, but when combining the algorithms we should also take account the errors assigned in the other algorithms, which does not pose any additional problems. For example, when determining $w_{11}$ in $f_1$, we should remember that in $f_2$ we determine a $\delta$ such that $w_{23}=v_3(1+\delta)$ with $\delta=\mathcal{O}(\epsilon_{\mathrm{machine}})$; in $f_1$ we should then choose $$w_{11}=u_1(1+\epsilon_1)(1+\epsilon_2)(1+\epsilon_3)/(1+\delta)=u_1(1+\mathcal{O}(\epsilon_{\mathrm{machine}})).$$ In calculating $\tilde{f}=(f_1\oplus f_2)\oplus f_3$, we also get errors in the additions, but we just assign one $(1+\mathcal{O}(\epsilon_{\mathrm{machine}})$ factor to every element of our matrix; since $(1+\mathcal{O}(\epsilon_{\mathrm{machine}})(1+\mathcal{O}(\epsilon_{\mathrm{machine}})=(1+\mathcal{O}(\epsilon_{\mathrm{machine}})$, this does not yield a problem and we can construct a matrix $W=[w_1|w_2]\in\mathbb{R}^{3\times2}$, such that $\tilde{f}([u|v])=f(W):=w_1\times w_2$ and $\frac{\|[u|v]-W\|}{\|[u|v]\|}=\mathcal{O}(\epsilon_{\mathrm{machine}})$. (Here $f$ is the exact cross product mapping.)
This is not, however, the most straightforward algorithm: clearly more straightforward is the algorithm $$\hat{f}:[u|v]\mapsto\begin{bmatrix}\phantom{-}[\mathrm{fl}(u_2)\otimes\mathrm{fl}(v_3)]\ominus[\mathrm{fl}(u_3)\otimes\mathrm{fl}(v_2)]\\ -[\mathrm{fl}(u_1)\otimes\mathrm{fl}(v_3)]\oplus[\mathrm{fl}(u_3)\otimes\mathrm{fl}(v_1)]\\ \phantom{-}[\mathrm{fl}(u_1)\otimes\mathrm{fl}(v_2)]\ominus[\mathrm{fl}(u_2)\otimes\mathrm{fl}(v_1)]\end{bmatrix}.$$
However, I was not able to show that this algorithm is backward stable, since I wasn't able to distribute the errors in such a way that we can conclude that for some $W=[w_1|w_2]\in\mathbb{R}^{3\times2}$, $\hat{f}([u|v])=f(W):=w_1\times w_2$ and $\frac{\|[u|v]-W\|}{\|[u|v]\|}=\mathcal{O}(\epsilon_{\mathrm{machine}})$. (Here $f$ is the exact cross product mapping.)
My question is whether this algorithm $\hat{f}$ is backward stable. If so, do you have a proof, or a reference to one? If it's not, why so?