#copyright 1997 by A. Robertson

print(`Version of March 25, 1997`):
lprint(``):
print(`This is SCHUR, A Maple package to investigate `):
print(`the asymptotic behavoir of the minimum number of`):
print(`Schur Triples,`):
print(`which accompanies the paper`):
print(`"On the asymptotic behavior of Schur triples."`):
print(`Graham proved that asymptotic behavoir of the`):
print(`minimum number of Schur Triples`):
print(`<=n^2/18 and wondered if it can`):
print(`be improved`):
lprint(``):
print(`The most current version is available on WWW at:`):
print(` http://www.math.temple.edu/~aaron .`):
print(`Please report all bugs to: aaron@math.temple.edu .`):
print(`All bugs or other comments used will be acknowledged in future`):
print(`versions.`):
lprint(``):
print(`For general help, and a list of the available functions,`):
print(` type "ezra();". `):
lprint(``):

ezra:=proc()
if args=NULL then
print(`This is SCHUR, A Maple package to investigate `):
print(`the asymptotic behavoir of the minimum number of`):
print(`Schur Triples,`):
print(`which accompanies the paper`):
print(`"On the asymptotic behavior of Schur triples."`):
print(`Graham proved that asymptotic behavoir of the`):
print(` minimum number of Schur Triples`):
print(`<=n^2/18 and wondered if it can`):
print(`be improved`):
lprint(``):
print(`Contains procedures: LBeqns,mat,SimpMax,SimpMin,UBeqns `):
fi:

if nops([args])=1 and op(1,[args])=LBeqns then
print(`LBeqns(In,A,B,C,type) outputs a list of lower bounds of the functions `):
print(`alpha(x,y,z)+alpha(x',y',z'), where (x,y,z) and (x',y',z') are`):
print(`from (0,1)^3`):
print(`It inputs`):
print(`In=list of the lengths of intervals`):
print(`A=list of the number of 0's in the corresponding interval`):
print(`B=list of the number of 1's in the corresponding interval`):
print(`C=list of the number of 01's in the corresponding interval`):
print(`type=1,2,or 3 depending on the format outputted (see the program`):
fi:

if nops([args])=1 and op(1,[args])=UBeqns then
print(`UBeqns(In,A,B,C) uses the same inputs as LBeqns.`):
print(`It outputs a list of upper bounds of the functions`):
print(`alpha(x,y,z)+alpha(x',y',z')`);
fi:

if nops([args])=1 and op(1,[args])=mat then
print(`mat(r) inputs the number of colors, r, and outputs`):
print(`the coefficient matrix.  For example mat(2) will give`):
print(`the matrix A in the paper.`):
fi:

if nops([args])=1 and (op(1,[args])=SimpMin or op(1,args)=SimpMax) then
print(`SimpMin(In,A,B,C) inputs numerical values`):
print(`for In,A,B, and C (cannot be symbolic)`);
print(`and uses the simplex method to output the minimum number of`);
print(`Schur Triples (alpha(0,0,0)+alpha(1,1,1)) based on the constraints`);
print(`computed using LBeqns and UBeqns.`);
print(`SimpMax(In,A,B,C) is the same as SimpMin except that it returns`);
print(`the maximum number of Schur Triples based on the constraints.`);
fi:

end:

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

LBeqns:=proc(In,A,B,C,type)
local i,j,X,Z,l,L,a,b,c,d,n:
L[0]:=0:

for i from 1 to nops(In) do
	a[i]:=A[i]:
	b[i]:=B[i]:
	c[i]:=C[i]:
	l[i]:=In[i]:
	L[i]:=sum(l[j],j=1..i):
	d[i]:=binomial(l[i],2)-binomial(a[i],2)-binomial(b[i],2)-c[i]:
od:

for i from 1 to 12 do
X[i]:=0:
od:

n:=L[nops(In)]:

for j from 1 to nops(In) do
	for i from 1 to j-1 do
		if L[i]+L[j] <= n then
			X[1]:=X[1]+a[i]*a[j]:
			X[2]:=X[2]+b[i]*b[j]:
			X[3]:=X[3]+a[i]*b[j]:
			X[4]:=X[4]+b[i]*a[j]:
		fi:

		if 2*L[i-1]+1 >= L[j] then
			X[5]:=X[5]+a[i]*a[j]:
			X[6]:=X[6]+b[i]*b[j]:
			X[7]:=X[7]+a[i]*b[j]:
			X[8]:=X[8]+b[i]*a[j]:
		fi:

		if 2*L[i] <= L[j-1] then
			X[9]:=X[9]+a[i]*a[j]:
			X[10]:=X[10]+b[i]*b[j]:
			X[11]:=X[11]+a[i]*b[j]:
			X[12]:=X[12]+b[i]*a[j]:
		fi:
	od:
od:

for i from 1 to nops(In) do
	if 2*L[i]-1 <=n then
		X[1]:=X[1]+binomial(a[i],2):
		X[2]:=X[2]+binomial(b[i],2):
		X[3]:=X[3]+c[i]:
		X[4]:=X[4]+d[i]:
	fi:

	if l[i]<=L[i-1] then
		X[5]:=X[5]+binomial(a[i],2):
		X[6]:=X[6]+binomial(b[i],2):
		X[7]:=X[7]+c[i]:
		X[8]:=X[8]+d[i]:
	fi:
od:

Z[1]:=alpha[0,0,0] + alpha[0,0,1]:
Z[2]:=alpha[1,1,0] + alpha[1,1,1]:
Z[3]:=alpha[0,1,0] + alpha[0,1,1]:
Z[4]:=alpha[1,0,0] + alpha[1,0,1]:
Z[5]:=alpha[0,0,0] + alpha[1,0,0]:
Z[6]:=alpha[0,1,1] + alpha[1,1,1]:
Z[7]:=alpha[0,0,1] + alpha[1,0,1]:
Z[8]:=alpha[0,1,0] + alpha[1,1,0]:
Z[9]:=alpha[0,0,0] + alpha[0,1,0]:
Z[10]:=alpha[1,0,1] + alpha[1,1,1]:
Z[11]:=alpha[0,0,1] + alpha[0,1,1]:
Z[12]:=alpha[1,0,0] + alpha[1,1,0]:

if type=1 then
{X[1]<=Z[1],X[2]<=Z[2],X[3]<=Z[3],X[4]<=Z[4],X[5]<=Z[5],X[6]<=Z[6],
X[7]<=Z[7],X[8]<=Z[8],X[9]<=Z[9],X[10]<=Z[10],X[11]<=Z[11],X[12]<=Z[12]}:
elif type=2 then
[X[1],X[2],X[3],X[4],X[5],X[6],X[7],X[8],X[9],X[10],X[11],X[12]]:
elif type=3 then
{X[1]=Z[1],X[2]=Z[2],X[3]=Z[3],X[4]=Z[4],X[5]=Z[5],X[6]=Z[6],
X[7]=Z[7],X[8]=Z[8],X[9]=Z[9],X[10]=Z[10],X[11]=Z[11],X[12]=Z[12]}:
fi:
end:

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

UBeqns:=proc(In,A,B,C)
local i,j,Y,Z,l,L,a,b,c,d,n,X:
L[0]:=0:
Y:=LBeqns2(In,A,B,C):

for i from 1 to nops(In) do
	a[i]:=A[i]:
	b[i]:=B[i]:
	c[i]:=C[i]:
	l[i]:=In[i]:
	L[i]:=sum(l[j],j=1..i):
	d[i]:=binomial(l[i],2)-binomial(a[i],2)-binomial(b[i],2)-c[i]:
od:

for i from 1 to 12 do
X[i]:=Y[i]:
od:

n:=L[nops(In)]:

for j from 1 to nops(In) do
	for i from 1 to j-1 do
		if L[i]+L[j] > n and L[i-1]+L[j-1]+2 <=n then
			X[1]:=X[1]+a[i]*a[j]:
			X[2]:=X[2]+b[i]*b[j]:
			X[3]:=X[3]+a[i]*b[j]:
			X[4]:=X[4]+b[i]*a[j]:
		fi:

		if 2*L[i-1]+1 < L[j]  and 2*L[i] >= L[j-1]+2 then
			X[5]:=X[5]+a[i]*a[j]:
			X[6]:=X[6]+b[i]*b[j]:
			X[7]:=X[7]+a[i]*b[j]:
			X[8]:=X[8]+b[i]*a[j]:
		fi:

		if 2*L[i] > L[j-1] and L[j]-2 >= 2*L[i-1]  then
			X[9]:=X[9]+a[i]*a[j]:
			X[10]:=X[10]+b[i]*b[j]:
			X[11]:=X[11]+a[i]*b[j]:
			X[12]:=X[12]+b[i]*a[j]:
		fi:
	od:
od:

for i from 1 to nops(In) do
	if 2*L[i-1] <= L[i]-2 then
		X[9]:=X[9]+binomial(a[i],2):
		X[10]:=X[10]+binomial(b[i],2):
		X[11]:=X[11]+c[i]:
		X[12]:=X[12]+d[i]:
	fi:

	if l[i]>L[i-1] then
		X[5]:=X[5]+binomial(a[i],2):
		X[6]:=X[6]+binomial(b[i],2):
		X[7]:=X[7]+c[i]:
		X[8]:=X[8]+d[i]:
	fi:
od:
Z[1]:=alpha[0,0,0] + alpha[0,0,1]:
Z[2]:=alpha[1,1,0] + alpha[1,1,1]:
Z[3]:=alpha[0,1,0] + alpha[0,1,1]:
Z[4]:=alpha[1,0,0] + alpha[1,0,1]:
Z[5]:=alpha[0,0,0] + alpha[1,0,0]:
Z[6]:=alpha[0,1,1] + alpha[1,1,1]:
Z[7]:=alpha[0,0,1] + alpha[1,0,1]:
Z[8]:=alpha[0,1,0] + alpha[1,1,0]:
Z[9]:=alpha[0,0,0] + alpha[0,1,0]:
Z[10]:=alpha[1,0,1] + alpha[1,1,1]:
Z[11]:=alpha[0,0,1] + alpha[0,1,1]:
Z[12]:=alpha[1,0,0] + alpha[1,1,0]:

{Z[1]<=X[1],Z[2]<=X[2],Z[3]<=X[3],Z[4]<=X[4],Z[5]<=X[5],Z[6]<=X[6],
Z[7]<=X[7],Z[8]<=X[8],Z[9]<=X[9],Z[10]<=X[10],Z[11]<=X[11],Z[12]<=X[12]}:

end:

##################################################3

SimpMin:=proc(In,A,B,C)
local f,set1,set2,constr,i,su,S,t:
with(simplex):
su:=0:

f:=alpha[0,0,0]+alpha[1,1,1]:
set1:=LBeqns(In,A,B,C):
set2:=UBeqns(In,A,B,C):
constr:=set1 union set2:

S:=minimize(f,constr,NONNEGATIVE):

for i from 1 to 8 do
  t:=op(i,S):

  if member(op(1,t),{alpha[0,0,0],alpha[1,1,1]})=true then
    su:=su+ op(2,t):
  fi:
od:
su:
end:

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

SimpMax:=proc(In,A,B,C)
local f,set1,set2,constr,i,su,S,t:
with(simplex):
su:=0:

f:=alpha[0,0,0]+alpha[1,1,1]:
set1:=LBeqns(In,A,B,C):
set2:=UBeqns(In,A,B,C):
constr:=set1 union set2:

S:=maximize(f,constr,NONNEGATIVE):

for i from 1 to 8 do
  t:=op(i,S):

  if member(op(1,t),{alpha[0,0,0],alpha[1,1,1]})=true then
    su:=su + op(2,t):
  fi:
od:
su:
end:

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

#given r colors outputs the matrix corresponding to the alpha's
mat:=proc(r)
local M,a,b,c,d,try1,try2,i,j,k,n,m,p,q,t,g,sq,L,count,w,ind,S,z:
n:=r^3-1:    # n+1 is the number of variables
p:=(r^2)*3:  #3 choices for which 2 places are same, r^2 for what is the same
M:=array(1..p,1..n+1):
for i from 1 to n+1 do
a[i]:=convert(i-1,base,r):
if nops(a[i]) < 3 then
  m[i]:=3-nops(a[i]):
  sq[i]:=seq(0,j=1..m[i]):
  b[i]:=[op(a[i]),sq[i]]:
else
  b[i]:=a[i]:
fi:
g[i]:=[seq(b[i][4-q],q=1..3)]:#re-ordering so they now are increasing
od:



L:=[[1,2],[1,3],[2,3]]:  #list of which entries can be the same
w:=p/3:                  #number of rows with each of the above the same
t:=1:

for z from 1 to p do

if z=1 or z=w+1 or z=2*w + 1 then
  S:={}:
fi:


  while t<>0 do
    for i from 1 to n+1 do
     if member(g[i],S)=false then
      try1:=g[i]:
      t:=0:
     fi:
    od:
  od:
  t:=1:

  for j from 1 to n+1 do
  try2:=g[j]:
  
    if z <= w then 
      ind:=L[1]:
      if try1[ind[1]]=try2[ind[1]] and try1[ind[2]]=try2[ind[2]] then
        M[z,j]:=1:
        S:=S union {try1,try2}:
      else
        M[z,j]:=0:
      fi:
      
    elif z > w and z <= 2*w then
      ind:=L[2]:
      if try1[ind[1]]=try2[ind[1]] and try1[ind[2]]=try2[ind[2]] then 
        M[z,j]:=1:
        S:=S union {try1,try2}:
      else
        M[z,j]:=0:
      fi:
         
    elif z > 2*w then
      ind:=L[3]:
      if try1[ind[1]]=try2[ind[1]] and try1[ind[2]]=try2[ind[2]] then 
        M[z,j]:=1:
        S:=S union {try1,try2}
      else
        M[z,j]:=0: 
      fi:      
    fi:

  od:
od:
M:
end:

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

bis:=proc(r)   #gives the b vector notation in Ax=b
local n,i,E,ans:
n:=r^2*3:
ans:=[]:
for i from 1 to n do
  ans:=[op(ans),E[i]]:
od:
end:
    
#######################################################