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