function search = Six_Round_Random_Offset(eta, number_searches)

%% Searches for 6-Round Bit-Commitment Based Coin-Flipping Protocols 
%% with bias less than 0.2499 for messages of dimension 2 
%% with the mesh offset by a randomly chosen parameter delta
%%
%% eta: precision 
%% number_searches: the number of searches ran with a new mesh offset 
%%
%% Worst case and average case tables are outputted displaying 
%% the efficiency of the strategies in the protocol filter.
%%
%% Questions? Contact Jamie at jamiesikora@gmail.com

d = eta; 

tic

Output = [];

minG1 = 7;
minG2 = 7;
minG3 = 7;
minG4 = 7;
minG5 = 7;
minG6 = 7;
minG7 = 7;
minG8 = 7;
minG9 = 7;
minG10 = 7;
minG11 = 7;
minG12 = 7;
minSDPA0 = 7;
minSDPA1 = 7;
minSDPB0 = 7;
minSDPB1 = 7;

totG1 = 0;
totG2 = 0;
totG3 = 0;
totG4 = 0;
totG5 = 0;
totG6 = 0;
totG7 = 0;
totG8 = 0;
totG9 = 0;
totG10 = 0;
totG11 = 0;
totG12 = 0;
totSDPA0 = 0;
totSDPA1 = 0;
totSDPB0 = 0;
totSDPB1 = 0;

for nom = 1:number_searches

nom

delta = rand(1)/100;

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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for c5 =  delta:prec:1
for c6 =  delta:prec:1-c5
for c8 =  delta:prec:1-c5-c6
c9 =  1-c5-c6-c8;
					   
beta0 = [c5,c6,c8,c9]'
iterbeta0 = iterbeta0 + 1;

beta_max = max(beta0);
					   
for d5 =  delta:prec:min(1, beta_max)
for d6 =  delta:prec:min(1-d5, beta_max)       
for d8 =  delta:prec:min(1-d5-d6, beta_max) 
d9 = 1-d5-d6-d8;
if d9 < beta_max + 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 =  delta:prec:1
for b6 =  delta:prec:1-b5
for b8 =  delta:prec:1-b5-b6
b9 =  1-b5-b6-b8;
					   
alpha0 = [b5,b6,b8,b9]';
alpha_max = max(alpha0);

for a5 =  delta:prec:min(1, alpha_max)
for a6 =  delta:prec:min(1-a5, alpha_max)       
for a8 =  delta:prec:min(1-a5-a6, alpha_max) 
a9 = 1-a5-a6-a8;
if a9 < alpha_max + 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;

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

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 = sedumi(A_conc, b_conc, -c_conc, K_conc, pars);
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 = sedumi(A_conc, b_conc, -c_conc, K_conc, pars);
		 G10 = x'*c_conc; % concave hull
if G10 < Best_Bias
iterG10 = iterG10 + 1;

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

A4 = kron(alpha0',eye(4)); % [A4 -eye(9)] gives constraint 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 = sedumi(A_Bob, b_Bob,-c0_Bob, K_Bob, pars);
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 = sedumi(A_Alice, b_Alice, -c_Alice, K_Alice, pars);
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 = sedumi(A_Bob, b_Bob,-c1_Bob, K_Bob, pars);
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 = sedumi(A_Alice, b_Alice, -c_Alice, K_Alice, pars);
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 % alpha0 generation
		 
end
end
end
end % alpha1 generation
		 
end
end
end % beta0 generation 

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

AS = iterbeta0beta1^2;

pG1 = 1-iterG1/AS;
pG2 = 1-iterG2/iterG1;
pG3 = 1-iterG3/iterG2;
pG4 = 1-iterG4/iterG3;
pG5 = 1-iterG5/iterG4;
pG6 = 1-iterG6/iterG5;
pG7 = 1-iterG7/iterG6;
pG8 = 1-iterG8/iterG7;
pG9 = 1-iterG9/iterG8; 
pG10 = 1-iterG10/iterG9; 
pSDPB0 = 1-iterSDPB0/iterG10; 
pG11 = 1-iterG11/iterSDPB0;
pSDPA0 = 1-iterSDPA0/iterG11; 
pSDPB1 = 1-iterSDPB1/iterSDPA0; 
pG12 = 1-iterG12/iterSDPB1; 
pSDPA1 = 1-iterSDPA1/iterG12; 
					
minG1 = min(minG1, pG1);					
minG2 = min(minG2, pG2);					
minG3 = min(minG3, pG3);
minG4 = min(minG4, pG4);					
minG5 = min(minG5, pG5);					
minG6 = min(minG6, pG6);					
minG7 = min(minG7, pG7);					
minG8 = min(minG8, pG8);					
minG9 = min(minG9, pG9);					
minG10 = min(minG10, pG10);					
minSDPB0 = min(minSDPB0, pSDPB0);					
minG11 = min(minG11, pG11);					
minSDPA0 = min(minSDPA0, pSDPA0);					
minSDPB1 = min(minSDPB1, pSDPB1);					
minG12 = min(minG12, pG12);					
minSDPA1 = min(minSDPA1, pSDPA1);

totG1 = totG1 + pG1;
totG2 = totG2 + pG2;
totG3 = totG3 + pG3;
totG4 = totG4 + pG4;
totG5 = totG5 + pG5;
totG6 = totG6 + pG6;
totG7 = totG7 + pG7;
totG8 = totG8 + pG8;
totG9 = totG9 + pG9;
totG10 = totG10 + pG10;
totG11 = totG11 + pG11;
totG12 = totG12 + pG12;
totSDPB0 = totSDPB0 + pSDPB0;
totSDPB1 = totSDPB1 + pSDPB1;
totSDPA0 = totSDPA0 + pSDPA0;
totSDPA1 = totSDPA1 + pSDPA1;

Total_Time = toc;

end
	
worst = [ ...
 minG1; ...			
 minG2; ...					
 minG3; ...
 minG4; ...					
 minG5; ...					
 minG6; ...					
 minG7; ...				
 minG8; ...					
 minG9; ...					
 minG10; ...	
 minSDPB0; ...					
 minG11; ...					
 minSDPA0; ...					
 minSDPB1; ...					
 minG12; ...				
 minSDPA1 ...
]

average = (1/number_searches)*[ ...   
 totG1; ...			
 totG2; ...					
 totG3; ...
 totG4; ...					
 totG5; ...					
 totG6; ...					
 totG7; ...				
 totG8; ...					
 totG9; ...					
 totG10; ...	
 totSDPB0; ...					
 totG11; ...					
 totSDPA0; ...					
 totSDPB1; ...					
 totG12; ...				
 totSDPA1 ...
]
					
return