%Tests kak1.m

global_defs;

example = 18
switch example
	case (1)
		U = rand_unitary(4,false);
	case (2)%U= Ax1
		U = kron(rand_unitary(2,true), id2);
	case (3)%U= 1xB
		U = kron(id2, rand_unitary(2,true));
	case (4)%U= AxB
		U = kron(rand_unitary(2,true), rand_unitary(2,true));
	case (5)%CNOT, sigx(0)^n(1) 
		U = [id2,ze2;
		ze2,sigx];
	case (6)%Controlled-U, A(0)^n(1)
		U = [id2,ze2;
		ze2,rand_unitary(2,true)];
	case (7)%Controlled-U, A(0)^nbar(1)
		U = [rand_unitary(2,true), ze2;
		ze2,id2];
	case (8)%multiplexor
		U = [rand_unitary(2,true), ze2;
		ze2,rand_unitary(2,true)];
	case (9)%swap, aka exchange operator 
		U =	ex_op;				
	case (10)%CNOT, sigx(1)^n(0)
		U =	[1,0,0,0;
			0,0,0,1;
			0,0,1,0;
			0,1,0,0];				
	case (11)
		cv = canonical_class_vec(randn(4,1))
		U = class_vec_to_mat(cv);
	case (12)
		cv1 = canonical_class_vec(randn(4,1))
		cv2 = canonical_class_vec(randn(4,1))
		U = class_vec_to_mat(cv1)*kron(
		rand_unitary(2,true), id2)*class_vec_to_mat(cv2);
	case (13)
		cv = canonical_class_vec(randn(4,1))
		U = [id2,ze2;
		ze2,rand_unitary(2,true)]*class_vec_to_mat(cv);
	case (14)
		a = randn(3,1);
		a = a/norm(a);
		cnot_a = [id2,ze2;ze2,sigx*a(1) + sigy*a(2) + sigz*a(3)];
		b = randn(3,1);
		%b = cross(a,b);%what if a and b are perpendicular?
		b = b/norm(b);		
		cnot_b = [id2,ze2;ze2,sigx*b(1) + sigy*b(2) + sigz*b(3)];
		%cv=[0;randn(2,1);0]; %cx=0
		%cv=[randn(2,1);0;0]; %cz=0
		cv=[randn();0; randn();0]; %cy=0
		U=cnot_a*class_vec_to_mat(cv)*cnot_b;
	case (15)
		a = randn(3,1);
		a = a/norm(a);
		cnot_a = cnot_a = [id2,ze2;ze2,sigx*a(1) + sigy*a(2) + sigz*a(3)];
		cv = canonical_class_vec(randn(4,1))
		U = cnot_a*class_vec_to_mat(cv);
	case (16)
		a = randn(3,1);
		a = a/norm(a);
		cnot_a = [id2,ze2;ze2,sigx*a(1) + sigy*a(2) + sigz*a(3)];
		cv=[randn(1,1);0;0;0];
		U=cnot_a*class_vec_to_mat(cv)*cnot_a;
	case (17)
		t1=randn();
		t2=randn();
		a=randn();
		b=randn();
		P1 = cos(t1)*id2 + i*sin(t1)*sigz;
		P2 = cos(t2)*id2 + i*sin(t2)*sigz;
		Pa = cos(a)*id2 + i*sin(a)*sigy;
		Pb = cos(b)*id2 + i*sin(b)*sigy;
		U=[id2,ze2;ze2,P1]*kron(Pa, Pb)*[id2,ze2;ze2,P2];
	case (18)
		a = [randn(); 0; randn()];
		b = [randn(); 0; randn()];
		alp = randn();
		bet = randn();
		mu = randn();
		nu = randn();
		%a(1)=0;
		%alp = 0;
		%alp = pi/2;
		%nu=0;
		mu=pi/4;
		a = a/norm(a); b = b/norm(b);
		cnot_a = [id2,ze2;ze2, sigx*a(1) + sigy*a(2) + sigz*a(3)];
		cnot_b = [id2,ze2;ze2, sigx*b(1) + sigy*b(2) + sigz*b(3)];
		rot_alp = kron(cos(alp)*id2 + i*sin(alp)*sigy, id2);
		rot_bet = kron(cos(bet)*id2 + i*sin(bet)*sigy, id2);
		e_mu = cos(mu) + i*sin(mu);
		d_mu = diag([e_mu, conj(e_mu), e_mu, conj(e_mu)]);
		e_nu = cos(nu) + i*sin(nu);
		d_nu = diag([e_nu, e_nu, conj(e_nu), conj(e_nu)]);
		d = d_nu*d_mu*[id2,ze2;ze2,i*sigz];
		U = cnot_a*rot_alp*d*rot_bet*cnot_b;
	otherwise
		error("example number is out of range");
end

[A1, A0, k, B1, B0]=kak1(U);

error=norm(kron(A1, A0)*
class_vec_to_mat(k)*kron(B1, B0)-U)

fi = fopen ("kak_out.txt", "w", "native");

fprintf(fi, "error = %10.4e\n", error);

fdisp(fi, "input U ="),fdisp(fi, U);

fdisp(fi, "A1 ="),fdisp(fi, A1);

fdisp(fi, "A0 ="),fdisp(fi, A0);

fdisp(fi, "B1 ="),fdisp(fi, B1);

fdisp(fi, "B0 ="),fdisp(fi, B0);

fdisp(fi, "class_vec ="),fdisp(fi, k);

fdisp(fi, "canonical_class_vec ="),fdisp(fi, canonical_class_vec(k));

fclose(fi);