4D rotation matrix to quaternions

Question

It is well-known that a 3D rotation can always be represented by a quaternion. It is less well-known that a 4D rotation can always be represented by two quaternions, sending a point \$p=(a,b,c,d)^T\$ represented as the quaternion \$a+bi+cj+dk\$ to the point \$p'=xpy\$ for an appropriate choice of unit quaternions \$x,y\$.

Let \$R\$ be a 4D rotation matrix. Then the corresponding \$x,y\$ can be obtained through a direct algorithm of Perez-Gracia and Thomas:

• Compute \$S=\frac{R-A_1RA_1-A_2RA_2-A_3RA_3}4\$ where $$A_1=\begin{bmatrix}0&0&0&-1\\0&0&-1&0\\0&1&0&0\\1&0&0&0\end{bmatrix}$$ $$A_2=\begin{bmatrix}0&0&1&0\\0&0&0&-1\\-1&0&0&0\\0&1&0&0\end{bmatrix}$$ $$A_3=\begin{bmatrix}0&-1&0&0\\1&0&0&0\\0&0&0&-1\\0&0&1&0\end{bmatrix}.$$
• Set \$Y=\frac S{(\det S)^{1/4}}\$ and \$X=RY^T\$. We have \$R=XY\$, where \$X\$ and \$Y\$ are the matrices of left- and right-isoclinic rotations respectively – very special matrices with at most four distinct magnitudes among their entries.
• Finally \$x\$ and \$y\$ are the left-most columns of \$X\$ and \$Y\$. The pair \$(-x,-y)\$ also corresponds to \$R\$, analogous to the 3D case.

\$S\$ above may have zero determinant, in which case one of the following alternative expressions should be used. These of course give different \$S\$ matrices, but the steps afterwards are unchanged: $$S=-\frac14(RA_1+A_1R+A_3RA_2-A_2RA_3)$$ $$S=-\frac14(RA_2+A_2R+A_1RA_3-A_3RA_1)$$ $$S=-\frac14(RA_3+A_3R+A_2RA_1-A_1RA_2)$$

Task

Given \$R\$ – an orthogonal (to floating-point error) 4×4 real matrix with determinant 1 – output a corresponding \$(x,y)\$ pair to within \$10^{-6}\$ entrywise relative error. You may use any reasonable and consistent formats for input and output, and any correct formula for the calculation (including the one in the previous section).

This is code-golf; fewest bytes wins.

Test cases

These are in the format \$R\to x\ y\$. The random cases were generated (and you can generate your own random cases) with this program.

[[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]] -> [1.0, 0.0, 0.0, 0.0] [1.0, 0.0, 0.0, 0.0]
[[-1, 0, 0, 0], [0, -1, 0, 0], [0, 0, -1, 0], [0, 0, 0, -1]] -> [1.0, 0.0, 0.0, 0.0] [-1.0, 0.0, 0.0, 0.0]
[[0, 0, 1, 0], [0, 0, 0, 1], [0, 1, 0, 0], [1, 0, 0, 0]] -> [1.0, 0.0, 0.0, 0.0] [0.0, 0.0, -1.0, 0.0]
[[0, 0, 1, 0], [-1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 0, -1]] -> [0.5, 0.5, 0.5, -0.5] [-0.5, -0.5, -0.5, -0.5]
[[0, 0, 0, -1], [1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0]] -> [0, -0.7071067811865476, 0, -0.7071067811865476] [-0.7071067811865476, 0, 0.7071067811865476, 0]
[[-0.6609058926555349, -0.30938449301078536, -0.6536493902436745, -0.2005669742797215], [-0.19737145221236216, 0.6931746532606626, 0.08266116424964398, -0.6882735952028355], [0.13703192282789176, -0.6502921800262064, 0.3686714537074865, -0.6499412988447983], [0.7109642487712476, 0.030169788832262925, -0.6557369638659809, -0.25224721015459695]] -> [0.560009784749661, 0.400037446835586, 0.272063553490942, 0.6725626401873549] [0.06637964559158858, 0.052592253308694104, 0.6737046111530368, 0.7341320688092741]
[[0.6544007065477859, -0.5813716570041381, -0.411235410397809, -0.2542678684713569], [0.6193101080207767, 0.7694074810291105, -0.005562224168405557, -0.1563207596643309], [0.38183012824424395, -0.12606111825246455, 0.23940297340390976, 0.8837423628837136], [-0.2059622745505508, 0.23265364607854452, -0.8795116345009956, 0.36043183434001774]] -> [0.9455738853169363, -0.26286929698266226, 0.18881691258838904, -0.033733271870118764] [0.5350303732855687, 0.7836340999984829, 0.31251919145079116, 0.044630168680125246]

---

This question and its 3D version were inspired by my answer to a MathsSE question. I needed to wrestle with 4D rotations when rendering images like these for this other MathsSE question which I have asked, and for which I have no answers:

Answer 1

Itr, 194 bytes

$R(0ää1~,0ä1~0,0 1áä,1 0ää)$A(0ä1á,0ää1~,1~0ää,0 1áä)$B(0 1~0ä,1 0ää,0ää1~,0ä1á)$C(RARA**-BRB**-CRC**-4÷RA*AR*+CRB**+BRC**-4÷~RB*BR*+ARC**+CRA**-4÷~RC*CR*+BRA**+ARB**-4÷~)äa¿×èäarr÷$Y(YRT*YT)µ0@

with newlines for readability:

$R
(0ää1~,0ä1~0,0 1áä,1 0ää)$A
(0ä1á,0ää1~,1~0ää,0 1áä)$B
(0 1~0ä,1 0ää,0ää1~,0ä1á)$C
(
 RARA**-BRB**-CRC**-4÷
 RA*AR*+CRB**+BRC**-4÷~
 RB*BR*+ARC**+CRA**-4÷~
 RC*CR*+BRA**+ARB**-4÷~
)
äa¿×è
äarr÷$Y
(YRT*YT)µ0@

input as matrix, output as pair of lists

For some test-cases \$(-x,-y)\$ is returned instead

online interpreter

version taking nested arrays as input

Explanation

Straight forward implementation of the algorithm in the question:

®                           ; convert (implicit) input to matrix (makes it easier to use the test cases)
$R                          ; store input in R
(0ää1~,0ä1~0,0 1áä,1 0ää)$A ; matrix A1 (0 1 -> literals, ä -> dup, ~ -> negate, á -> over)
(0ä1á,0ää1~,1~0ää,0 1áä)$B  ; matrix A2
(0 1~0ä,1 0ää,0ää1~,0ä1á)$C ; matrix A3
(                           ; array containing the three possibilities for S
 RARA**-BRB**-CRC**-4÷
 RA*AR*+CRB**+BRC**-4÷~
 RB*BR*+ARC**+CRA**-4÷~
 RC*CR*+BRA**+ARB**-4÷~
)
                            ; get all elements of the array with non-zero determinant:
äa                          ; vector of determinants of the elements
  ¿                         ; is nonzero (vectorized)
   ×                        ; repeat elements of array that many times, picks elements satisfying condition
è                           ; get the first element of that array "S"
äarr÷$Y                     ; compute Y ( S/ sqrt(sqrt(det( S )))) 
(YRT*YT)                    ; pair containing X^T=Y*R^T and Y^T
µ0@                         ; replace elements with their first row (get first column of X and Y)