# restart; read `c:/Users/Asus/Google Drive/Aek/NewSchur/Elaine.txt`;

with(Optimization):

######################################
# Section1: Optimization
# Proposition 20 Section 3
######################################

# Try: A11A(R1,R2,B,B1);
A11A := proc(R1,R2,B,B1) option remember;
local g1,g,A;

A := (B-g)/2*g1;
expand(subs({g1=B1,g=R1-R2},A));
end:

# Try: Delta11A();
Delta11A := proc() option remember;
local R,R1,R2,B,B1,C,Delta;

Delta := 1/2*(R*B/2+R*B1+B*R1+2*A11A(R1,R2,B,B1));
C := {0 <=R1, R1 <= R, R <=1/2}
union {B1<=R1, R1>=R2, B1<=(R1-R2)/2};

Delta := subs({R2=R-R1,B=1-R,B1=1/2-R1},Delta);
C := subs({R2=R-R1,B=1-R,B1=1/2-R1},C);

Maximize(Delta,C);
end:


# Try: A11B(R1,R2,B,B1);
A11B := proc(R1,R2,B,B1) option remember;
local g1,g,A;

A := (B-g)/2*g1;
expand(subs({g1=R1,g=R1-R2},A));
end:

# Try: Delta11B();
Delta11B := proc() option remember;
local R,R1,R2,B,B1,C,Delta;

Delta := 1/2*(R*B/2+R*B1+B*R1+2*A11B(R1,R2,B,B1));
C := {0 <=R1, R1 <= R, R <=1/2}
union {B1>=R1, R1 >=R2, R1<=(R1-R2)/2};

Delta := subs({R2=R-R1, B=1-R, B1=1/2-R1},Delta);
C := subs({R2=R-R1, B=1-R, B1=1/2-R1},C);

Maximize(Delta,C);
end:


# Try: A11C(R1,R2,B,B1);
A11C := proc(R1,R2,B,B1) option remember;
local g1,g,A;

A := (B-g)/2*g1;
expand(subs({g1=R1,g=R2-R1},A));
end:

# Try: Delta11C();
Delta11C := proc() option remember;
local R,R1,R2,B,B1,C,Delta;

Delta := 1/2*(R*B/2+R*B1+B*R1+2*A11C(R1,R2,B,B1));
C := {0 <=R1, R1 <= R, R <=1/2}
union {B1>=R1, R1 <=R2, R1<=(R2-R1)/2};

Delta := subs({R2=R-R1, B=1-R, B1=1/2-R1},Delta);
C := subs({R2=R-R1, B=1-R, B1=1/2-R1},C);

Maximize(Delta,C);
end:


# Try: A12(R1,R2,B,B1);
A12 := proc(R1,R2,B,B1) option remember;
local g1,g,A;

A := (B-g)/2*g1;
expand(subs({g1=B/4,g=B/2},A));
end:


# Try: Delta12A();
Delta12A := proc() option remember;
local R,R1,R2,B,B1,C,Delta;

Delta := 1/2*(R*B/2+R*B1+B*R1+2*A12(R1,R2,B,B1));
C := {0 <=R1, R1 <= R, R <=1/2}
union {B1<=R1, R1>=R2, B1 >=(R1-R2)/2};

Delta := subs({R2=R-R1,B=1-R,B1=1/2-R1},Delta);
C := subs({R2=R-R1,B=1-R,B1=1/2-R1},C);

Maximize(Delta,C);
end:


# Try: Delta12B();
Delta12B := proc() option remember;
local R,R1,R2,B,B1,C,Delta;

Delta := 1/2*(R*B/2+R*B1+B*R1+2*A12(R1,R2,B,B1));
C := {0 <=R1, R1 <= R, R <=1/2}
union {B1>=R1, R1 >=R2, R1>=(R1-R2)/2};

Delta := subs({R2=R-R1, B=1-R, B1=1/2-R1},Delta);
C := subs({R2=R-R1, B=1-R, B1=1/2-R1},C);

Maximize(Delta,C);
end:


# Try: Delta12C();
Delta12C := proc() option remember;
local R,R1,R2,B,B1,C,Delta;

Delta := 1/2*(R*B/2+R*B1+B*R1+2*A12(R1,R2,B,B1));
C := {0 <=R1, R1 <= R, R <=1/2}
union {B1>=R1, R1<=R2, R1>=(R2-R1)/2};

Delta := subs({R2=R-R1, B=1-R, B1=1/2-R1},Delta);
C := subs({R2=R-R1, B=1-R, B1=1/2-R1},C);

Maximize(Delta,C);
end:

#####################################################

# Try: A21A(R1,R2,B,B1);
A21A := proc(R1,R2,B,B1) option remember;
local g1,g,A;

A := 1/2*(B-g1)*g1-g^2/8;
expand(subs({g1=B1,g=R1-R2},A));
end:

# Try: Delta21A();
Delta21A := proc() option remember;
local R,R1,R2,B,B1,C,Delta;

Delta := 1/2*(R*B/2+R*B1+B*R1+2*A21A(R1,R2,B,B1));
C := {0 <=R1, R1 <= R, R <=1/2}
union {B1<=R1, R1>=R2, (R1-R2)/2 <= B1, B1<=R/2};

Delta := subs({R2=R-R1,B=1-R,B1=1/2-R1},Delta);
C := subs({R2=R-R1,B=1-R,B1=1/2-R1},C);

Maximize(Delta,C);
end:



# Try: A21B(R1,R2,B,B1);
A21B := proc(R1,R2,B,B1) option remember;
local g1,g,A;

A := 1/2*(B-g1)*g1-g^2/8;
expand(subs({g1=R1,g=R1-R2},A));
end:

# Try: Delta21B();
Delta21B := proc() option remember;
local R,R1,R2,B,B1,C,Delta;

Delta := 1/2*(R*B/2+R*B1+B*R1+2*A21B(R1,R2,B,B1));
C := {0 <=R1, R1 <= R, R <=1/2}
union {B1>=R1, R1>=R2, (R1-R2)/2 <= R1, R1<=R/2};

Delta := subs({R2=R-R1,B=1-R,B1=1/2-R1},Delta);
C := subs({R2=R-R1,B=1-R,B1=1/2-R1},C);

Maximize(Delta,C);
end:





# Try: A21C(R1,R2,B,B1);
A21C := proc(R1,R2,B,B1) option remember;
local g1,g,A;

A := 1/2*(B-g1)*g1-g^2/8;
expand(subs({g1=R1,g=R2-R1},A));
end:

# Try: Delta21C();
Delta21C := proc() option remember;
local R,R1,R2,B,B1,C,Delta;

Delta := 1/2*(R*B/2+R*B1+B*R1+2*A21C(R1,R2,B,B1));
C := {0 <=R1, R1 <= R, R <=1/2}
union {B1>=R1, R1<=R2, (R2-R1)/2 <= R1, R1<=R/2};

Delta := subs({R2=R-R1,B=1-R,B1=1/2-R1},Delta);
C := subs({R2=R-R1,B=1-R,B1=1/2-R1},C);

Maximize(Delta,C);
end:




# Try: A22A(R,R1,R2,B,B1);
A22A := proc(R,R1,R2,B,B1) option remember;
local g1,g,A;

A := 1/2*(B-g1)*g1-g^2/8;
expand(subs({g1=R/2,g=R1-R2},A));
end:

# Try: Delta22A();
Delta22A := proc() option remember;
local R,R1,R2,B,B1,C,Delta;

Delta := 1/2*(R*B/2+R*B1+B*R1+2*A22A(R,R1,R2,B,B1));
C := {0 <=R1, R1 <= R, R <=1/2}
union {B1<=R1, R1>=R2, (R1-R2)/2 <= R/2, R/2<=B1};

Delta := subs({R2=R-R1,B=1-R,B1=1/2-R1},Delta);
C := subs({R2=R-R1,B=1-R,B1=1/2-R1},C);

Maximize(Delta,C);
end:





# Try: A22B(R,R1,R2,B,B1);
A22B := proc(R,R1,R2,B,B1) option remember;
local g1,g,A;

A := 1/2*(B-g1)*g1-g^2/8;
expand(subs({g1=R/2,g=R1-R2},A));
end:

# Try: Delta22B();
Delta22B := proc() option remember;
local R,R1,R2,B,B1,C,Delta;

Delta := 1/2*(R*B/2+R*B1+B*R1+2*A22B(R,R1,R2,B,B1));
C := {0 <=R1, R1 <= R, R <=1/2}
union {B1>=R1, R1>=R2, (R1-R2)/2 <= R/2, R/2<=R1};

Delta := subs({R2=R-R1,B=1-R,B1=1/2-R1},Delta);
C := subs({R2=R-R1,B=1-R,B1=1/2-R1},C);

Maximize(Delta,C);
end:


# Try: A22C(R,R1,R2,B,B1);
A22C := proc(R,R1,R2,B,B1) option remember;
local g1,g,A;

A := 1/2*(B-g1)*g1-g^2/8;
expand(subs({g1=R/2,g=R2-R1},A));
end:

# Try: Delta22C();
Delta22C := proc() option remember;
local R,R1,R2,B,B1,C,Delta;

Delta := 1/2*(R*B/2+R*B1+B*R1+2*A22C(R,R1,R2,B,B1));
C := {0 <=R1, R1 <= R, R <=1/2}
union {B1>=R1, R1<=R2, (R2-R1)/2 <= R/2, R/2<=R1};

Delta := subs({R2=R-R1,B=1-R,B1=1/2-R1},Delta);
C := subs({R2=R-R1,B=1-R,B1=1/2-R1},C);

Maximize(Delta,C);
end:


################################################

# Try: A3A(R,R1,R2,B,B1);
A3A := proc(R,R1,R2,B,B1) option remember;
local g1,g,A;

A := (1/2-g1)*g1-g^2/8-R^2/8;
expand(subs({g1=B1,g=R1-R2},A));
end:

# Try: Delta3A();
Delta3A := proc() option remember;
local R,R1,R2,B,B1,C,Delta;

Delta := 1/2*(R*B/2+R*B1+B*R1+2*A3A(R,R1,R2,B,B1));
C := {0 <=R1, R1 <= R, R <=1/2}
union {B1<=R1, R1>=R2, R/2 <= B1};
 
Delta := subs({R2=R-R1,B=1-R,B1=1/2-R1},Delta);
C := subs({R2=R-R1,B=1-R,B1=1/2-R1},C);

Maximize(Delta,C);
end:





# Try: A3B(R,R1,R2,B,B1);
A3B := proc(R,R1,R2,B,B1) option remember;
local g1,g,A;

A := (1/2-g1)*g1-g^2/8-R^2/8;
expand(subs({g1=R1,g=R1-R2},A));
end:

# Try: Delta3B();
Delta3B := proc() option remember;
local R,R1,R2,B,B1,C,Delta;

Delta := 1/2*(R*B/2+R*B1+B*R1+2*A3B(R,R1,R2,B,B1));
C := {0 <=R1, R1 <= R, R <=1/2}
union {B1>=R1, R1>=R2, R/2 <= R1};

Delta := subs({R2=R-R1,B=1-R,B1=1/2-R1},Delta);
C := subs({R2=R-R1,B=1-R,B1=1/2-R1},C);

Maximize(Delta,C);
end:



# Try: A3C(R,R1,R2,B,B1);
A3C := proc(R,R1,R2,B,B1) option remember;
local g1,g,A;

A := (1/2-g1)*g1-g^2/8-R^2/8;
expand(subs({g1=R1,g=R2-R1},A));
end:

# Try: Delta3C();
Delta3C := proc() option remember;
local R,R1,R2,B,B1,C,Delta;

Delta := 1/2*(R*B/2+R*B1+B*R1+2*A3C(R,R1,R2,B,B1));
C := {0 <=R1, R1 <= R, R <=1/2}
union {B1>=R1, R1<=R2, R/2 <= R1};

Delta := subs({R2=R-R1,B=1-R,B1=1/2-R1},Delta);
C := subs({R2=R-R1,B=1-R,B1=1/2-R1},C);

Maximize(Delta,C);
end:







##################################
# Section 2: Optimization x+ay=z,
# Proposition 28, Section 4
##################################
# Input: a 
# Output: The number of maximum number of 
# non-mono-chromatic triples, ratio of R1 and R. 
# Try: Winner(13); GenDelta3A(13);
Winner := proc(a) option remember;
[evalf(1/(2*a)*(1-1/(a^2+2*a+3))),
evalf((a+1)/(a^2+2*a+3)),evalf((a+2)/(a^2+2*a+3))];
end:


# General $a$ of case3A
# Try: GenA3A(a,R,R1,R2,B,B1);
GenA3A := proc(a,R,R1,R2,B,B1) option remember;
local g1,g,A;

A := a/2*(1/a-g1)*g1-g^2/(4*a)-R^2/(4*a);
expand(subs({g1=B1,g=R1-R2},A));
end:

# Try: GenDelta3A(3);
GenDelta3A := proc(a) option remember;
local R,R1,R2,B,B1,C,Delta;

Delta := 1/2*(R*B/a+R*B1+B*R1+2*GenA3A(a,R,R1,R2,B,B1));
C := {0 <=R1, R1 <= R, R <=1/2}
union {B1<=R1, R1>=R2, R/a <= B1};
 
Delta := subs({R2=R-R1,B=1-R,B1=1/a-R1},Delta);
C := subs({R2=R-R1,B=1-R,B1=1/a-R1},C);

Maximize(Delta,C);
end:




###############################
# Section 3: Conjectures of the 
# problems in section 5
###############################

# Input: a,b and n the length of the interval
# Output: the interval which give the 
# minimum mono-chormatic triples a*x+b*y=a*z.
# Try: POne(1,2,20);
POne := proc(a,b,n) option remember;
if gcd(a,b) > 1 then
	ERROR("a,b");
fi:	
Find1([1],a,b,n,[n^2,{}]);
end:

# Try: Find1([1],1,1,5,[5^2,{}]);
Find1 := proc(C,a,b,n,Ans) option remember;
local T1;

if nops(C) = n then
	if Trip(C,a,b) < Ans[1] then
		return([Trip(C,a,b),{C}]);
	elif Trip(C,a,b) = Ans[1] then
		return([Ans[1], Ans[2] union {C}]);
	else
		return(Ans);	
	fi:	
elif Trip(C,a,b) > Ans[1] then
	return(Ans);	
fi:

T1 := Find1([op(C),1],a,b,n,Ans);  
Find1([op(C),2],a,b,n,T1);  

end:

# Count number of mono-chromatic triples 
# a*x+b*y=a*z on interval C 
# seq(Trip([1$k],1,2),k=1..12);
Trip := proc(C,a,b) option remember;
local x,n,count;

n := nops(C);
if n=0 then
	return(0);
fi:
	
count := Trip( [op(1..n-1,C)],a,b);

for x from 1 to n-1 do 
	if C[x] = C[n] and a*(n-x) mod b = 0
	and a*(n-x)/b <=n and C[a*(n-x)/b] = C[n] then
		count := count+1;
	fi:	
od:

count;
end:


######################################################
# Input: n the length of the interval
# Output: the interval which give the 
# minimum mono-chormatic triples x+y+w=z.
# Try: PTwo(20);
PTwo := proc(n) option remember;
Find2([1],n,[n^2,{}]);
end:

# Try: Find2([1],5,[5^2,{}]);
Find2 := proc(C,n,Ans) option remember;
local T1;

if Quad(C) > Ans[1] then
	return(Ans);
elif nops(C) = n then
	if Quad(C) < Ans[1] then
		return([Quad(C),{C}]);
	elif Quad(C) = Ans[1] then
		return([Ans[1], Ans[2] union {C}]);
	else
		return(Ans);	
	fi:	
fi:

T1 := Find2([op(C),1],n,Ans);  
Find2([op(C),2],n,T1);  

end:

# Count number of mono-chromatic triples 
# x+y+w=z on interval C 
# seq(Quad([1$k]),k=1..12);
Quad := proc(C) option remember;
local x,y,n,count;

n := nops(C);
if n=0 then
	return(0);
fi:
	
count := Quad( [op(1..n-1,C)]);

for x from 1 to n-1 do 
	for y from 1 to n-x-1 do
		if nops({C[x],C[y],C[n-x-y],C[n]})=1  then
			count := count+1;
		fi:	
	od:	
od:

count;
end:





########################################################
# Input: n the length of the interval
# Output: the interval which give the maximum 
# rainbow triples x+y=z under 3-coloring.
# Try: PThree(13);
# [seq(PThree(k)[1],k=1..14)];
PThree := proc(n) option remember;
Find3([1],n,[0,{}]);
end:

# Try: Find3([1],5,[0,{}]);
Find3 := proc(C,n,Ans) option remember;
local T1;

if nops(C) = n then
	if Rainbow(C) > Ans[1] then
		return([Rainbow(C),{C}]);
	elif Rainbow(C) = Ans[1] then
		return([Ans[1], Ans[2] union {C}]);
	else
		return(Ans);	
	fi:	
fi:

T1 := Find3([op(C),1],n,Ans);  
T1 := Find3([op(C),2],n,T1);  
Find3([op(C),3],n,T1);
end:

# Count number of rainbow triples 
# x+y=z on interval C 
# seq(Rainbow([1$k]),k=1..12);
Rainbow := proc(C) option remember;
local x,n,count;

n := nops(C);
if n=0 then
	return(0);
fi:
	
count := Rainbow( [op(1..n-1,C)]);

for x from 1 to n/2 do 
	if {C[x],C[n-x],C[n]} ={1,2,3} then
		count := count+1;
	fi:	
od:

count;
end:



#####################################
# Just in case, 
# but I did not see the pattern
#####################################
# Now with equation x+y=w+z
# Try: Find22([1],5,[5^2,{}]);
Find22 := proc(C,n,Ans) option remember;
local T1;

if TwoPairs(C) > Ans[1] then
	return(Ans);
elif nops(C) = n then
	if TwoPairs(C) < Ans[1] then
		return([TwoPairs(C),{C}]);
	elif TwoPairs(C) = Ans[1] then
		return([Ans[1], Ans[2] union {C}]);
	else
		return(Ans);	
	fi:	
fi:

T1 := Find22([op(C),1],n,Ans);  
Find22([op(C),2],n,T1);  

end:

# Count number of mono-chromatic triples 
# x+y=w+z on interval C 
# seq(TwoPairs([1$k]),k=1..12);
TwoPairs := proc(C) option remember;
local x,y,n,count;

n := nops(C);
if n=0 then
	return(0);
fi:
	
count := TwoPairs( [op(1..n-1,C)]);

for x from 1 to n do 
	for y from n+1-x to n do
		if nops({C[x],C[y],C[x+y-n],C[n]})=1  then
			count := count+1;
		fi:	
	od:	
od:

count;
end: