%Check explicit expressions for LO-RHS invariant
%of 2-bit circuit with 3 dressed cnots, 
%in its diagonalized form,
%assuming b1 perp span(a1, c1)
%and {c0,b0,a0} is a dependent set

clear;

parametrize=1;

global_defs;
for j=1:3

%use randn()instead of randn() to get negatives
a0= get_unit_vec(randn(3,1));
b0= get_unit_vec(randn(3,1));
c0= get_unit_vec(randn()*a0 +randn()*b0);	
	
a1= get_unit_vec(randn(3,1));
b1= get_unit_vec(cross(a1,randn(3,1)));
c1= get_unit_vec(cross(b1,randn(3,1)));

if(parametrize)
	if(-(cross(a0,b0)'*cross(b0,c0))*(b0'*c0)<0)
		printf("a0 changed\n");
		a0=-a0;
	endif
endif


cross_ab1= cross(a1,b1);
cross_bc1 = cross(b1,c1);
%cross_abc1= cross(cross_ab1,c1);
%cross_abcc1 = cross(cross_abc1,c1);
%cross_abb1 = cross(cross_ab1,b1);
%cross_abbc1 = cross(cross_abb1,c1);
%cross_abbcc1 = cross(cross_abbc1,c1);
%cross_bcc1 = cross(cross_bc1,c1);

cross_ab0= cross(a0,b0);
cross_bc0 = cross(b0,c0);
%cross_abc0= cross(cross_ab0,c0);
%cross_abcc0 = cross(cross_abc0,c0);
cross_abb0 = cross(cross_ab0,b0);
cross_abbc0 = cross(cross_abb0,c0);
cross_abbcc0 = cross(cross_abbc0,c0);
cross_bcc0 = cross(cross_bc0,c0);

%vol_0 = cross_ab0'*c0
vol_1 = cross_ab1'*c1;

ac1 = a1'*c1;
ab0 = a0'*b0;
bc0 = b0'*c0;

sin_bet2 = norm(cross_bc0);
eta = norm(cross_abbc0);
xxo = ac1*(cross_abb0'*c0);
yyo = ab0*bc0*vol_1;


mu = [-ac1*eta, -ab0*ac1*sin_bet2, (cross_abb0'*c0)*vol_1];

nu = [ab0*vol_1*sin_bet2, vol_1*eta, ab0*bc0*ac1];

if(parametrize)
	cos1=ab0;
	sin1=norm(cross_ab0);
	cos2=bc0;
	sin2=norm(cross_bc0);
	xi=sign(cross_ab0'*cross_bc0);
	xi2=sign(cos2);
	cos_bet = ac1;
	sin_bet = vol_1;
	alt_mu = [-cos_bet*sin1*cos2*xi2,-cos_bet*cos1*sin2,xi*sin_bet*sin1*sin2];
		err_mu = norm(mu-alt_mu)
	alt_nu = [sin_bet*cos1*sin2,sin_bet*sin1*cos2*xi2,cos_bet*cos1*cos2];
		err_nu = norm(nu-alt_nu)
	alt_xxo = cos_bet*sin1*sin2*xi;
		err_xxo = abs(xxo-alt_xxo)
	alt_yyo = sin_bet*cos1*cos2;
		err_yyo = abs(yyo-alt_yyo)
endif

gx0= c0;
gy0= cross_bc0/sin_bet2;
gz0= -cross_bcc0/sin_bet2;

gx1= c1;
gy1= cross_bc1;
gz1= b1;

v1 = gy0;
v2 = -cross_abbcc0/eta;
v3 = gx0;
if(parametrize)
	alt_v2 = gz0;
	err_v2 = norm(v2-alt_v2)
endif

u1 = cross_abbc0/eta;
u2 = gz0;
u3 = gx0;
if(parametrize)
	alt_u1 = gy0;
	err_u1 = norm(u1-alt_u1)
endif

inv1 = dr11(c1,c0)*dr11(b1,b0)*dr11(a1,a0)*dr11(-a1,-a0)*dr11(-b1,-b0)*dr11(-c1,-c0);
inv2 = (-xxo -i*yyo)*eye(4) + nu(1)*kron(sig(gx1),sig(v1)) + nu(2)*kron(sig(gy1),sig(v2)) + nu(3)*kron(sig(gz1),sig(v3)) +i*(mu(1)*kron(sig(gx1),sig(u1)) + mu(2)*kron(sig(gy1),sig(u2))+ mu(3)*kron(sig(gz1),sig(u3)));
err_3cnot = norm(inv1 - inv2)

endfor