function	[U1,U0] = factor_SU2pow2_matrix(U)

%This function checks whether a matrix U is a
%tensor product of two SU(2) matrices,
%U=U1 otimes U0.
%If it is, the function finds U1 and U0.


%Start by assuming that 
%U =
%[x*U0, y*U0]
%[-y'*U0, x'*U0]
%Will check validity of this assumption at the end

x_sq = U(3:4,3)'*U(1:2,1);
y_sq = -U(3:4,1)'*U(1:2,3);
yh_x = U(1:2,3)'*U(1:2,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

U1=[x,y;-y',x'];
if(abs(x)>abs(y))
	U0 = U(1:2,1:2)/x;
else
	U0 = U(1:2,3:4)/y;
end


%the moment of truth
kk = norm(U-kron(U1,U0));


if(kk>1e-9)
	error("U doesnt factor into tensor product");
end