%Check explicit expressions for LO-RHS invariant
%of 2-bit circuit with 1,2 and 3 dressed cnots.
%ckt_invar = circuit invariant

clear;

global_defs;

%use randn()instead of rand() to get negatives

a0= get_unit_vec(randn(3,1));
a1= get_unit_vec(randn(3,1));
b0= get_unit_vec(randn(3,1));
b1= get_unit_vec(randn(3,1));
c0= get_unit_vec(randn(3,1));
c1= get_unit_vec(randn(3,1));

cross_ab0 = cross(a0,b0);
cross_ab1 = cross(a1,b1);
cross_bc0 = cross(b0,c0);
cross_bc1 = cross(b1,c1);
cross_bcc0 = cross(cross_bc0,c0);
cross_bcc1 = cross(cross_bc1,c1);
cross_abb0 = cross(cross_ab0,b0);
cross_abb1 = cross(cross_ab1,b1);
cross_abc0 = cross(cross_ab0,c0);
cross_abc1 = cross(cross_ab1,c1);
cross_abcc0 = cross(cross_abc0,c0);
cross_abcc1 = cross(cross_abc1,c1);
cross_abbc0 = cross(cross_abb0,c0);
cross_abbc1 = cross(cross_abb1,c1);
cross_abbcc0 = cross(cross_abbc0,c0);
cross_abbcc1 = cross(cross_abbc1,c1);

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

error_1cnot = norm(
-kron(sig(a1),sig(a0))
-dr11(a1,a0)*dr11(-a1,-a0))

error_2cnots = norm(
(a1'*b1)*(a0'*b0)*eye(4)
-kron(sig(cross_abb1), sig(cross_abb0))
+ i*(a0'*b0)*kron(sig(cross_ab1), sig(b0))
+ i*(a1'*b1)*kron(sig(b1), sig(cross_ab0))
-dr11(b1,b0)*dr11(a1,a0)*dr11(-a1,-a0)*dr11(-b1,-b0))


lam3r = (cross_abb1'*c1)*(cross_abb0'*c0);
lam3i = -((a0'*b0)*(b0'*c0)*vol_1 +(a1'*b1)*(b1'*c1)*vol_0);

L3r = (-(a1'*b1)*(a0'*b0)*kron(sig(c1),sig(c0))
+(a0'*b0)*(b0'*c0)*kron(sig(cross_abc1),sig(c0))
+(a1'*b1)*(b1'*c1)*kron(sig(c1),sig(cross_abc0))
+(a1'*b1)*vol_0*kron(sig(cross_bc1),sig(c0))
+(a0'*b0)*vol_1*kron(sig(c1),sig(cross_bc0))
-kron(sig(cross_abbcc1),sig(cross_abbcc0)));

L3i = ((a0'*b0)*kron(sig(cross_abcc1),sig(cross_bcc0))
+(a1'*b1)*kron(sig(cross_bcc1),sig(cross_abcc0))
+(cross_abb0'*c0)*kron(sig(cross_abbc1),sig(c0))
+(cross_abb1'*c1)*kron(sig(c1),sig(cross_abbc0)));

 
error_3cnots = norm(eye(4)*(lam3r + i*lam3i) + L3r + i*L3i
-dr11(c1,c0)*dr11(b1,b0)*dr11(a1,a0)*dr11(-a1,-a0)*dr11(-b1,-b0)*dr11(-c1,-c0))

L3rG = Gamma_rep(L3r);
L3iG = Gamma_rep(L3i);

err_commutator=norm(L3r*L3i-L3i*L3r)
err_right_commutator = norm(L3rG*transpose(L3iG)-L3iG*transpose(L3rG))
err_left_commutator = norm(transpose(L3rG)*L3iG-transpose(L3iG)*L3rG)