function [A,B] = SO4_to_SU2xSU2(Q)

%Given a matrix Q\in SO(4), this function
%finds matrices A, B\in SU(2) such that Q = M'*kron(A, B)*M

global_declarations;

if(norm(Q*transpose(Q)-eye(4))>1e-8)
	error("input matrix for SO4_to_SU2xSU2 is not orthogonal");
end

if(abs(det(Q)-1)>1e-8)
	error("input matrix for SO4_to_SU2xSU2 is not special");
end

U = uni2ort*Q*uni2ort';

%Start by assuming that 
%U =
%[a11*B, a12*B]
%[-a12'*B, a11'*B]
%Will check validity of this assumption at the end

A = zeros(2,2);
B = zeros(2,2);

a11_sq = (U(1:2,1:2)*U(3:4,3:4)')(1,1);
a12_sq = -(U(1:2,3:4)*U(3:4,1:2)')(1,1);
a11_a12H = (U(1:2,1:2)*U(1:2,3:4)')(1,1);
A(1,1) = sqrt(a11_sq);
A(1,2) = sqrt(a12_sq);

%sqrt ambiguous in its sign, 
%so the values just
%assigned to A(1,1) and
%to A(1,2) may have the wrong sign.
%Fix relative sign:
if(abs(A(1,1)*A(1,2)'-a11_a12H)>1e-9)
	A(1,2)=-A(1,2);
end

A(2,1) = -A(1,2)';
A(2,2) = A(1,1)';
if(abs(A(1,1))>abs(A(1,2)))
	B = U(1:2,1:2)/A(1,1);
else
	B = U(1:2,3:4)/A(1,2); 
end

%the moment of truth
if(norm(kron(A,B)-U)>1e-9)
	%U
	%kron(A,B)
	error("the error in  SO4_to_SU2xSU2 is too large")
end