function [U2,U1,U0] = factor_SU2pow3_matrix(U)
%This function checks whether a matrix U is a
%tensor product of three SU(2) matrices,
%U=U2 otimes U1 otimes U0,
%If it is, the function finds U2, U1 and U0.
%Start by assuming that
%U =
%[x*U1U0, y*U1U0]
%[-y'*U1U0, x'*U1U0]
%where U1U0 = U1 otimes U0.
%Will check validity of this assumption at the end.
x_sq = U(5:8,5)'*U(1:4,1);
y_sq = -U(5:8,1)'*U(1:4,5);
yh_x = U(1:4,5)'*U(1:4,1);
x = sqrt(x_sq);
y = sqrt(y_sq);
%sqrt ambiguous in its sign,
%so the values just
%assigned to x and
%to y may have the wrong sign.
%Fix relative sign:
if(abs(y'*x-yh_x)>1e-9)
x=-x;
end
U2=[x,y;-y',x'];
if(abs(x)>abs(y))
U1U0 = U(1:4,1:4)/x;
else
U1U0 = U(1:4,5:8)/y;
end
[U1,U0] = factor_SU2pow2_matrix(U1U0);
%the moment of truth
kk = norm(U-kron(U2,kron(U1,U0)));
if(kk>1e-9)
error("U doesnt factor into tensor product");
end