function search = Six_Round_Search_d2(eta_min, eta_max)

%% Searches for 6-Round Bit-Commitment Based Coin-Flipping Protocols 
%% with bias less than 0.2499 for messages of dimension 2
%%
%% Input: (eta_min, eta_max), where 
%% eta_min is the coarsest mesh precision, e.g., 2
%% eta_max is the finest mesh precision, e.g., 5
%%
%% Each column of Output is the result of a search 
%% (The first entry is the precision, the second is how many 
%% made it past the symmetry reductions, and then the rest are how
%% many protocols made it past each strategy in the filter.)
%% 
%% Note: the program also outputs values of beta0 to see if progress is being made.
%%
%% Questions? Contact Jamie at jamiesikora@gmail.com

max_gap = 0;

tic

Output = [];

for d = eta_min:eta_max

E = 1;

b = 0.75-0.0001;

MD = kron(eye(2),[1 1]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Search Parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

prec = 1/d;

format('long')

Best_Bias = b;					% Searching for protocols with bias less than this

iter = 0;

pars.fid = 0;

e = [1 1];

iterbeta0 = 0;
iterbeta0beta1 = 0;

iterG1 = 0;
iterG2 = 0;
iterG3 = 0;
iterG4 = 0;
iterG5 = 0;
iterG6 = 0;
iterG7 = 0;
iterG8 = 0;
iterG9 = 0;
iterG10 = 0;
iterG11 = 0;
iterG12 = 0;
iterSDPA0 = 0;
iterSDPA1 = 0;
iterSDPB0 = 0;
iterSDPB1 = 0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Alice SDP (Protocol Independent)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

D = [];
for i = 1:4
temp = sparse(4,4);
temp(i,i) = 1;
D = [D; vec(temp)'];
end

e= [1 1];

K_Alice.l = 2 + 8; K_Alice.s = [4 4 4 4  4 4 4 4];

A_Alice = [ ...
e sparse(1,8) sparse(1,8*16); ... % tr(c) = 1
-kron(ones(2,1),eye(2)) kron(eye(2), kron(eye(2),e)) sparse(4,8*16); ... % tr_A2(s^y1) = c
sparse(16,2) kron(eye(2), kron([1;1], eye(4))) -kron(eye(4), [D D]) ...  % sum_a diag(sigma^(a,y)) = s^y1
];

b_Alice = [1; sparse(16+4,1)];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Alice concave hull (protocol independent)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

K_conc.l = 4; K_conc.s = [4 4 4 4  4 4 4 4];

A_conc = [ ...
eye(4) sparse(4,8*16); ... % setting the p^y1 = v
kron(ones(4,1), eye(4)) -kron(eye(4), [D D]) ... % sum_a diag(sigma^(a,y)) = s^y1
];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Bob SDP (Protocol Independent)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

K_Bob.l = 4+16+8;
K_Bob.s = [4 4];

A2 = kron(eye(2),e); % sum of the q_{x_1}'s is 1.

A3 = [kron(eye(2),kron(e',eye(2))) -kron(eye(4),kron(eye(2),e))];
% tr_B2(p_x1x2) = q_x1 constraint

b_Bob = [1;1;sparse(8+4+4+4+4,1)];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Protocol Generation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Can permute A (and B) such that A1 has:
% First block is non-decreasing
% First element in each block is non-decreasing.

for c5 =  0:prec:1
for c6 =  0:prec:1-c5
for c8 =  0:prec:1-c5-c6
c9 =  1-c5-c6-c8;
if c9 > max([c5,c6,c8]) - 0.0001
					   
beta0 = [c5,c6,c8,c9]'
max_gap;
					   
for d5 =  0:prec:min(1,c9)
for d6 =  0:prec:min(1-d5,c9)       
for d8 =  0:prec:min(1-d5-d6,c9) 
d9 = 1-d5-d6-d8;
if d9 < c9 + 0.0001
					   
beta1 = [d5,d6,d8,d9]';
					   
iterbeta0beta1 = iterbeta0beta1 + 1;

G1 = 0.5 + 0.5*sqrt(beta0)'*sqrt(beta1);
if G1 < Best_Bias

for b5 =  0:prec:1
for b6 =  0:prec:1-b5
for b8 =  0:prec:1-b5-b6
b9 =  1-b5-b6-b8;
if b9 > max([b5,b6,b8]) - 0.0001
					   
alpha0 = [b5,b6,b8,b9]';
					   
for a5 =  0:prec:min(1,b9)
for a6 =  0:prec:min(1-a5,b9)       
for a8 =  0:prec:min(1-a5-a6,b9) 
a9 = 1-a5-a6-a8;
if a9 < b9 + 0.0001
					   
alpha1 = [a5,a6,a8,a9]';

iterG1 = iterG1 + 1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%% 	Filter (Known Cheating Strategies) 	%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

rho0 = MD*alpha0;
rho1 = MD*alpha1;					   

G2 = 0.5 + 0.25 * norm(rho0 - rho1,1);
if G2 < Best_Bias
iterG2 = iterG2 + 1;
					   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

sigma0 = MD*beta0;
sigma1 = MD*beta1;					   

w=0;
v=0;		  
for counter = 1:4		  
if alpha1(counter) > alpha0(counter)
v = v + alpha1(counter); 
else w = w + alpha0(counter);
end
end	
	
G3 = 0.25*(w+v+sqrt((w-v)^2 +4*w*v*(sqrt(sigma0)'*sqrt(sigma1))^2));
if G3 < Best_Bias
iterG3 = iterG3 + 1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

G4 = (0.5 + 0.5*sqrt(rho0)'*sqrt(rho1)) * (0.5 + 0.25 * norm(sigma0 - sigma1,1));
if G4 < Best_Bias
iterG4 = iterG4 + 1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
					
m=0;
n=0;		  
for counter = 1:4		  
if beta1(counter) > beta0(counter)
n = n + beta1(counter); 
else m = m + beta0(counter);
end
end	
G5 = 0.25*(m+n+sqrt((m-n)^2 +4*m*n*(sqrt(alpha0)'*sqrt(alpha1))^2));
if G5 < Best_Bias
iterG5 = iterG5 + 1;
					
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
			
m=0;
n=0;		  
for counter = 1:2		  
if sigma1(counter) > sigma0(counter)
     n = n + sigma1(counter); 
else m = m + sigma0(counter);
end
end	
	  
G6 = 0.25*(m+n+sqrt((m-n)^2 +4*m*n*(sqrt(rho0)'*sqrt(rho1))^2));
if G6 < Best_Bias
iterG6 = iterG6 + 1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

ATemp = 0.5*w*sqrt(sigma0)*sqrt(sigma0)' + 0.5*v*sqrt(sigma1)*sqrt(sigma1)';
n = null(ATemp - G3*eye(2)); % change G3 possibly if name changes
S = size(n);
if S(2) > 0
  c = n(:,1);
  else
  [c, e1] = eigs(ATemp,1);
end
c = abs(c); 

for ic = 1:2
  if sigma0(ic) + sigma1(ic) < prec/2
    c(ic) = 0;
  end
end

c = c/norm(c,2); c = c.*c;

p2 = sparse(4,4);
for counter1 = 1:4 % x
for counter2 = 1:4 % y    
if counter2 < 2.5
    y1 = 1;
else 
    y1 = 2;
end

    if alpha1(counter1) > alpha0(counter1)
        % g = 1
      if sigma1(y1) > 0.0001
          p2(counter1, counter2) = c(y1)*beta1(counter2)/sigma1(y1);
      else
          p2(counter1, counter2) = c(y1)/2;
      end
    else
        % g = 0
      if sigma0(y1) > 0.0001
          p2(counter1, counter2) = c(y1)*beta0(counter2)/sigma0(y1);
      else
          p2(counter1, counter2) = c(y1)/2;
      end
    end
end
end 

p0 = 0;
p1 = 0;
for i = 1:4
 p0 = p0 + alpha0(i)*p2(i,:)';
 p1 = p1 + alpha1(i)*p2(i,:)';
end

G7 = 0.5*((sqrt(p0)'*sqrt(beta0))^2+(sqrt(p1)'*sqrt(beta1))^2);
if G7 < Best_Bias
iterG7 = iterG7 + 1;

%[rand; Output]

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

ATemp = 0.5*w*sqrt(sigma1)*sqrt(sigma1)' + 0.5*v*sqrt(sigma0)*sqrt(sigma0)';
n = null(ATemp - G3*eye(2)); % change G3 possibly if name changes
S = size(n);
if S(2) > 0
  c = n(:,1);
  else
  [c, e1] = eigs(ATemp,1);
end
c = abs(c); 

for ic = 1:2
  if sigma0(ic) + sigma1(ic) < prec/2
    c(ic) = 0;
  end
end
c = c/norm(c,2); c = c.*c;

p2 = sparse(4,4);
for counter1 = 1:4 % x
for counter2 = 1:4 % y    
if counter2 < 2.5
    y1 = 1;
else 
    y1 = 2;
end

    if alpha1(counter1) > alpha0(counter1)
        % g = 1
      if sigma0(y1) > 0.0001
          p2(counter1, counter2) = c(y1)*beta0(counter2)/sigma0(y1);
      else
          p2(counter1, counter2) = c(y1)/2;
      end
    else
        % g = 0
      if sigma1(y1) > 0.0001
          p2(counter1, counter2) = c(y1)*beta1(counter2)/sigma1(y1);
      else
          p2(counter1, counter2) = c(y1)/2;
      end
    end
end
end 

p0 = 0;
p1 = 0;
for i = 1:4
 p0 = p0 + alpha0(i)*p2(i,:)';
 p1 = p1 + alpha1(i)*p2(i,:)';
end

G8 = 0.5*((sqrt(p0)'*sqrt(beta1))^2+(sqrt(p1)'*sqrt(beta0))^2);
if G8 < Best_Bias
iterG8 = iterG8 + 1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

	  ATemp = sparse(4,4);
	  for counter = 1:4
	  if beta1(counter) > beta0(counter)
        ATemp = ATemp + 0.5 * beta1(counter) * sqrt(alpha1)*sqrt(alpha1)';
        else ATemp = ATemp + 0.5 * beta0(counter) * sqrt(alpha0)*sqrt(alpha0)';
	  end
	  end	  

	  
m0 = vec(sqrt(alpha0)*sqrt(alpha0)');
m1 = vec(sqrt(alpha1)*sqrt(alpha1)');

c_conc = 0.5*[sparse(4,1); ...
beta0(1)*m0; beta1(1)*m1; ...
beta0(2)*m0; beta1(2)*m1; ...
beta0(3)*m0; beta1(3)*m1; ...
beta0(4)*m0; beta1(4)*m1 ...
];

n = null(ATemp - G5*eye(4)); % Change this if there are numerical errors
S = size(n);
if S(2) > 0
  v1 = n(:,1);
  else
	  [v1, e1] = eigs(ATemp,1);
end
v1 = abs(v1); 
for ic = 1:4
  if alpha0(ic) + alpha1(ic) < prec/2
    v1(ic) = 0;
  end
end

v1 = v1/norm(v1,2); v1 = v1.*v1;
b_conc = [v1; sparse(16,1)];

[x,y,INFO] = sedumi(A_conc, b_conc, -c_conc, K_conc, pars);
gap = abs(x'*c_conc + y'*b_conc);
if gap > max_gap
    max_gap = gap
end  
if max([INFO.pinf, INFO.dinf]) > 0.5
   [INFO.pinf, INFO.dinf]					
   error('SDP Wrong')
end 
G9 = x'*c_conc; % concave hull
if G9 < Best_Bias
iterG9 = iterG9 + 1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	 
c_conc = 0.5*[sparse(4,1); ...
					   beta1(1)*m0; beta0(1)*m1; ...
					   beta1(2)*m0; beta0(2)*m1; ...
					   beta1(3)*m0; beta0(3)*m1; ...
					   beta1(4)*m0; beta0(4)*m1 ...
					   ];
					
		 ATemp = sparse(4,4);
		 for counter = 1:4
		 if beta1(counter) > beta0(counter)
		 ATemp = ATemp + 0.5 * beta1(counter) * sqrt(alpha0)*sqrt(alpha0)';
		 else ATemp = ATemp + 0.5 * beta0(counter) * sqrt(alpha1)*sqrt(alpha1)';
		 end
		 end
					
n = null(ATemp - G5*eye(4)); % Change this if there are numerical errors
		 S = size(n);
		 if S(2) > 0
		 v1 = n(:,1);
		 else
						   [v1, e1] = eigs(ATemp,1);
		 end
		 v1 = abs(v1); 
         
for ic = 1:4
  if alpha0(ic) + alpha1(ic) < prec/2
    v1(ic) = 0;
  end
end
         v1 = v1/norm(v1,2); v1 = v1.*v1;
		 b_conc = [v1; sparse(16,1)];
		 
[x,y,INFO] = sedumi(A_conc, b_conc, -c_conc, K_conc, pars);
gap = abs(x'*c_conc + y'*b_conc);
if gap > max_gap
    max_gap = gap
end  
if max([INFO.pinf, INFO.dinf]) > 0.5
   [INFO.pinf, INFO.dinf]					
   error('SDP Wrong')
end 
		 G10 = x'*c_conc; % concave hull
if G10 < Best_Bias
iterG10 = iterG10 + 1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

A4 = kron(alpha0',eye(4)); % [A4 -eye(9)] gives constriant r_0 = sum_x alpha0_x p_x
A5 = kron(alpha1',eye(4));

A_Bob = [ ...
A2 sparse(2,16) sparse(2,8) sparse(2,2*16); ... % sum_y1 [q_x1]_y1 = 1
A3 sparse(8,8) sparse(8,2*16); ... % tr_B2(p_x1x2) = q_x1 constraint
sparse(4,4) A4 -eye(4) sparse(4,4) sparse(4,2*16); ... % r_0 = sum_x alpha0_x p_x
sparse(4,4) A5 sparse(4,4) -eye(4) sparse(4,2*16); ... % r_1 = sum_x alpha1_x p_x
sparse(4,4) sparse(4,16) -eye(4) sparse(4,4) D sparse(4,16); ...
sparse(4,4) sparse(4,16) sparse(4,4) -eye(4) sparse(4,16) D];

c0_Bob = 0.5*[sparse(4+16+8,1); vec(sqrt(beta0)*sqrt(beta0)'); vec(sqrt(beta1)*sqrt(beta1)')];

[x,y,INFO] = sedumi(A_Bob, b_Bob,-c0_Bob, K_Bob, pars);
gap = abs(x'*c0_Bob + y'*b_Bob);
if gap > max_gap
    max_gap = gap
end 
if max([INFO.pinf, INFO.dinf]) > 0.5
   [INFO.pinf, INFO.dinf]					
   error('SDP Wrong')
end 
SDPB0 = x'*c0_Bob;
if SDPB0 < Best_Bias
iterSDPB0 = iterSDPB0 + 1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

G11 = 1/(2*SDPB0);
if G11 < Best_Bias
iterG11 = iterG11 + 1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

c_Alice = 0.5*[sparse(2+8,1); ...
beta0(1)*m0; beta1(1)*m1; ...
beta0(2)*m0; beta1(2)*m1; ...
beta0(3)*m0; beta1(3)*m1; ...
beta0(4)*m0; beta1(4)*m1 ...
];

[x,y,INFO] = sedumi(A_Alice, b_Alice, -c_Alice, K_Alice, pars);
gap = abs(x'*c_Alice + y'*b_Alice);
if gap > max_gap
    max_gap = gap
end 
if max([INFO.pinf, INFO.dinf]) > 0.5
   [INFO.pinf, INFO.dinf]					
   error('SDP Wrong')
end 
temp = x'*c_Alice;
SDPA0 = temp(1,1);
if SDPA0 < Best_Bias
iterSDPA0 = iterSDPA0 + 1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

c1_Bob = 0.5*[sparse(4+16+8,1); vec(sqrt(beta1)*sqrt(beta1)'); vec(sqrt(beta0)*sqrt(beta0)')];

[x,y,INFO] = sedumi(A_Bob, b_Bob,-c1_Bob, K_Bob, pars);
gap = abs(x'*c1_Bob + y'*b_Bob);
if gap > max_gap
    max_gap = gap
end 
if max([INFO.pinf, INFO.dinf]) > 0.5
   [INFO.pinf, INFO.dinf] 					
   error('SDP Wrong')
end 
SDPB1 = x'*c1_Bob;
if SDPB1 < Best_Bias
iterSDPB1 = iterSDPB1 + 1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

G12 = 1/(2*SDPB1);
if G12 < Best_Bias
iterG12 = iterG12 + 1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

c_Alice = 0.5*[sparse(2+8,1); ...
beta1(1)*m0; beta0(1)*m1; ...
beta1(2)*m0; beta0(2)*m1; ...
beta1(3)*m0; beta0(3)*m1; ...
beta1(4)*m0; beta0(4)*m1 ...
];
		 
[x,y,INFO] = sedumi(A_Alice, b_Alice, -c_Alice, K_Alice, pars);
gap = abs(x'*c_Alice + y'*b_Alice);
if gap > max_gap
    max_gap = gap
end    
if max([INFO.pinf, INFO.dinf]) > 0.5
   [INFO.pinf, INFO.dinf]
	error('SDP Wrong')
end 
temp = x'*c_Alice;
SDPA1 = temp(1,1);
if SDPA1 < Best_Bias
iterSDPA1 = iterSDPA1 + 1;

alpha0
alpha1
beta0
beta1
error('New Protocol Found')

%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%	Filter Endifs
%%%%%%%%%%%%%%%%%%%%%%%%%

end % G1
end % G2
end % G3
end % G4
end % G5
end % G6
end % G7
end % G8
end % G9
end % G10
end % G11
end % G12

end % SDPA0
end % SDPA1
end % SDPB0
end % SDPB1

%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%	Protocol Generation Endifs
%%%%%%%%%%%%%%%%%%%%%%%%%

end
end
end
end % alpha0 generation
		 
end
end
end
end % alpha1 generation
		 
end
end
end
end % beta0 generation 

end
end
end % beta1 generation
end
		 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Output
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%TP = iterbeta0^4;
AS = iterbeta0beta1^2;

iter = [ ...
%TP; ...
d; ...
AS; ...
iterG1; ...
iterG2; ...
iterG3; ...
iterG4; ...
iterG5; ...
iterG6; ...
iterG7; ...
iterG8; ...
iterG9; ...
iterG10; ...
iterSDPB0; ...
iterG11; ...
iterSDPA0; ...
iterSDPB1; ...
iterG12; ...
iterSDPA1
];

Output = [Output, iter]

Total_Time = toc;

% max_gap % used for testing the SDP duality gap

end

return