%This subroutine checks Lemma 3, Eq.67 of
%version 2 of arXiv:0706.0479
%``How to Compile Some NAND Formula Evaluators",
%by R.R. Tucci
%The Lemma states an approximate
%matrix factorization that originates from
%the Campbell-Baker-Hausdorff (CBH) expansion.

%inputs:
%g = 1/sqrt(2)
g= .04;
mat_case = 1;

rt2 = sqrt(2);
for j=1:2 
	g = g + .01
	rand("seed", .1);
	
	switch(mat_case)
		case 1
			A = g*rand(4:4);
			A = A + A';
			B = g*rand(4:4);
			%B = B + B';
		case 2
			A = g*rt2*[0, 0, 0, 0;
			0, 0, 0, 1;
	 		0, 0, 0, 0;
			0, 1, 0, 0];
			
			B = g*rt2*[0, 0, 0, 0;
			0, 0, 0, 0;
	 		0, 0, 0, 0;
			0, 0, 0, 1];
	endswitch 
	
	F = expm(i*[A,B';B,zeros(4)]);	
		
	[V, D] = eig(A); %A = V*D*V';
	
	DVec = diag(D);
	
	sincD = zeros(4,4);
	expiD = zeros(4,4);
	for k=1:4
		sincD(k,k) = my_sinc(DVec(k)/2);
		expiD(k,k) = exp(i*DVec(k)/2);
	endfor
	
	L = B*V*sincD*expiD*V';
	BB = B'*B;
	AA = A*A;
	
	switch(mat_case)
		case 1
			F0 = [expm(i*A),zeros(4,4); zeros(4,4), eye(4)];
			F1 = expm([i/12*(BB*A + A*BB)+(1/24)*(- BB*AA + AA*BB), i*L'; i*L, -i*B*A*B'/6]    );
		case 2
		 	F0 = [expm(i*A + i/12*(BB*A + A*BB)),zeros(4,4); zeros(4,4), eye(4)];
		 	F1 = expm([zeros(4,4), i*L'; i*L, zeros(4,4)]    );
	endswitch
	
	error = norm(F - F0*F1, "fro")
endfor