function   [vecs_out, U1, U0] = deflate_dcnots(vecs_in)

% d0    t0    t0    a0         b0f  a0f  [U0]
% |     |     |     |    -->   |     |    
% d1    d1    a1    a1         b1f  a1f  [U1]
%A string of 2 controlled-Us can be easily expressed 
%as the LHS of the above identity. Thus, this function
%can be used to reduce two controlled-Us to 2 DC-CNOTs

d0=vecs_in(1:3,1);
%t0=vecs_in(1:3,2);
t0=vecs_in(1:3,3);
a0=vecs_in(1:3,4);

d1=vecs_in(4:6,1);
%d1=vecs_in(4:6,2);
%a1=vecs_in(4:6,3);
a1=vecs_in(4:6,4);

%cross_ad0 = cross(a0,d0);
cross_at0 = cross(a0,t0);
cross_td0 = cross(t0,d0);
cross_atd0 = cross(cross_at0,d0);
cross_att0 = cross(cross_at0,t0);
cross_atdd0 = cross(cross_atd0,d0);
cross_attd0 = cross(cross_att0,d0);

cross_ad1 = cross(a1,d1);
cross_add1 = cross(cross_ad1,d1);

at0 = a0'*t0;
td0 = t0'*d0;
ad1 = a1'*d1;

fy1 = cross_ad1;
sin_phi1 = norm(fy1);
fy1 = fy1/sin_phi1;

fy0 = d0;

c1_c0 = -at0*td0 + ad1*cross_att0'*d0;
s1_s0 = (cross_at0'*d0)*sin_phi1;

hx1 = cross_add1/sin_phi1;
hz1 = d1;

eta = norm(cross_atd0);
hx0 = cross_atd0/eta;
hz0 = cross_atdd0/eta;
cos_2 = cross_attd0'*hz0;
sin_2 = cross_attd0'*hx0;
cos_1 = cross_td0'*hz0;
sin_1 = cross_td0'*hx0;

M = [(-at0*sin_1 +ad1*sin_2),(-at0*cos_1 +ad1*cos_2);0,-sin_phi1*eta];

[s1_c0, c1_s0, fx0, fz0, fx1, fz1]= diag_ckt_invar2_aux(M, hx0, hz0, hx1, hz1);


c_sum = c1_c0 - s1_s0;
c_dif = c1_c0 + s1_s0;
s_sum = s1_c0 + c1_s0;
s_dif = s1_c0 - c1_s0;
alp_sum = atan2(s_sum,c_sum);
alp_dif = atan2(s_dif,c_dif);
alp1 = (alp_sum + alp_dif)/2;
alp0 = (alp_sum - alp_dif)/2;

b0f= fx0;
a0f = fx0*cos(alp0)- fy0*sin(alp0);

b1f= fx1;
a1f = fx1*cos(alp1)- fy1*sin(alp1);

vecs_out=[b0f,a0f;b1f,a1f];

%%%%%%%%%%%%%%%%%% now find U1 and U0
[U1,U0] = factor_SU2pow2_matrix(
dr11(a1f,a0f)*dr11(b1f,b0f)*dr11(d1,d0)*dr11(d1,t0)*dr11(a1,t0)*dr11(a1,a0));