#copyright 1998 by A. Robertson and D. Zeilberger

print(`Version of September 25, 1998`):
lprint(``):
print(`This is GENRON, A Maple package accompanying `):
print(`written by Aaron Robertson and Doron Zeilberger`):
print(` The lastest version of this package and of the article`):
print(`	are available from http://www.math.temple.edu/~zeilberg`):
print(`	and http://www.math.temple.edu/~aaron`):
lprint(``):
print(`Please report all bugs to: zeilberg@math.temple.edu .`):
print(`or 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();". For specific help type "ezra(procedure_name)" `):
lprint(``):
with(simplex): 
 
 
ezra:=proc()
if args=NULL then
print(`This is RON, A Maple package to investigate Ron Graham's`):
print(`problem about finding`):
print(`the number of monochromatic solns of x+ay=z`):
lprint(``):
print(`Contains the procedures:`):
print(`katan, schur, B, and Ron`):
fi:
  
if nops([args])=1 and op(1,[args])=katan then
print(`katan(n,m) finds the smallest possible number of`):
print(`monochromatic solutions of x+my=z`):
print(`that a 0-1 vector of length n can have`):
fi:
 
 
 
 
if nops([args])=1 and op(1,[args])=Ron then
print(`Ron(list,m) finds the asymptotics of schur(a,m)/n^2,`):
print(`when n goes to infinity`):
print(`where a is [0^(list[1]*n),1^((list[2]-list[1])*n), ...,]`):
print(`where list[n+1] is 1`):
print(`For example to find the limit of`):
print(` schur([0^{4/11*n},1^{6/11*n},0^{1/11*n}])/n^2,1),type`):
print( ` Ron([4,11,10/11],1) `):
fi:
 
 
if nops([args])=1 and op(1,[args])=B then
print(`B(list,m) finds the number of monochromatic`):
print(`solutions of x+my=z`):
print(`of the vector [0^list[1],1^list[2],0^list[3], ...]`):
fi:

 
if nops([args])=1 and op(1,[args])=schur then
print(`schur(a,m) finds the number of monochromatic`):
print(`solutions in the 0-1 vector a.`):
print(`In other words: the number of 1<=i<j<=n`):
print(` such that a[i]=a[j]=a[i+mj].`):
print(`The input a is given in terms of a list`):
print(`For example schur([0,0,0,1],1) should give 1`):
fi:
end:
 
###################################################
################################################### 
#looking for monochromatic solns (i,j,i+mj)
#i.e. monochromatics solns of x+my=z
#a is the list of the coloring
#i<j

schur:=proc(a,m)
local co,i,j,n:
 
co:=0:
 
n:=nops(a) :
 
for i from 1 to n do
for j from i+1 to (n-i)/m do
 if op(i,a)=op(j,a) and op(j,a)=op(i+m*j,a) then
  co:=co+1:
 fi:
od:
od: 
co: 
end:


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

#set of all 0-1 vectors of length n

vects:=proc(n)
local mu,gu,i:

if n=1 then
 RETURN({[1],[0]}):
fi:
 
if n=0 then
 RETURN({[]}):
fi:
 
if n<0 then
ERROR(`n>=0`):
fi:
 
mu:=vects(n-1):
gu:={}:
 
for i from 1 to nops(mu) do
gu:=gu union {[op(op(i,mu)),0],[op(op(i,mu)),1]}:
od:
gu:
end:

##################################################
 
katan:=proc(n,m)
local i,pu,ku,lu,gu,mu:
 
gu:=n^2:
 
mu:=vects(n):
 
lu:={}:
for i from 1 to nops(mu) do
ku:=op(i,mu):
pu:=schur(ku,m):
 
if pu=gu then
 gu:=pu:
 lu:=lu union {ku}:
fi:
 
 
if pu<gu then
 gu:=pu:
 lu:={ku}:
fi:
 
od:
 
gu,lu:
end:

##############################################
 
katanad:=proc(L,m)
local n,mu,gu:
gu:=[]:
for n from 1 to L do
print(`n=`,n):
mu:=katan(n,m):
print(mu):
gu:=[op(gu),mu]:
od:
gu:
end:
 
###############################################
 
B1:=proc(n,m)
local gu,i:
gu:=[seq(0,i=1..n)]:
schur(gu,m):
end:
 
###############################################
 
B:=proc(resh,m)
local gu,j,i,n:
n:=nops(resh):
gu:=[]:
for j from 1 to trunc(n/2) do
gu:=[op(gu),seq(0,i=1..op(2*j-1,resh) ) , seq(1,i=1..op(2*j,resh) ) ]:
od:
 
if not type(n/2, integer) then
gu:=[op(gu),seq(0,i=1..op(n,resh) ) ]:
fi:
schur(gu,m):
end:
 
##############################################
 
Br:=proc(resh,m):
evalf(B(resh,m)/convert(resh,`+`)^2):
end:

############################################## 
 
#f111(a,b,m) finds the number of triples (i,j,i+j) s.t. a<=i<j<i+mj<=b
f111:=proc(a,b,m):
 
if not (0<=a and a<b ) then
 ERROR(`0<=a<b`):
fi:
 
if b-(m+1)*a<=0 then
 RETURN(0):
else
 RETURN((b-(m+1)*a)^2 / (2*m*(m+1)) ):
fi:
end:
################################################# 
 
#f112(a,b,c,d,m) finds the number of triples (i,j,i+j) s.t. a<=i<j<=b
#and c<=i+mj<=d
 
f112:=proc(a,b,c,d,m)
local gu:
 
 
if not (0<=a and a<b and b<c and c<d) then
 ERROR(`0<=a<b<c<d`):
fi:
 
#Case I
 if c>=(m+1)*b then
  RETURN(0):
 fi:
 
#Case II
 if a+m*b<=c and c<(m+1)*b then
  gu:=((m+1)*b-c)^2 / (2*m*(m+1)):
  if (m+1)*b-d>=0 then
   gu:=gu-((m+1)*b-d)^2 / (2*m*(m+1)):
  fi:
  RETURN(gu):
 fi:
 
#case III
if (m+1)*a<=c and c<a+m*b then
#case III.1
  if d>=(m+1)*b then
    RETURN((b-a)^2/2-(c-(m+1)*a)^2 / (2*m*(m+1))):
  fi:
 
#case III.2
  if a+m*b<=d and d<(m+1)*b then
         RETURN((b-a)^2/2-(c-(m+1)*a)^2 / (2*m*(m+1)) - ((m+1)*b-d)^2 /
(2*m*(m+1))):
  fi:
   
#case III.3
  if (m+1)*a<=d and d<a+m*b then
         RETURN((d-(m+1)*a)^2 / (2*m*(m+1)) - (c-(m+1)*a)^2 / (2*m*(m+1))):
  fi:
   
fi:
 
#case IV
 
if c<(m+1)*a then
#case IV.1
  if d<=(m+1)*a then
     RETURN(0):
  fi:
#case IV.2
  if (m+1)*a<d and d<=a+m*b then
   RETURN((d-(m+1)*a)^2/(2*m*(m+1))):
  fi:
#case IV.3
  if a+m*b<d and d<=(m+1)*b then
    RETURN((b-a)^2/2-((m+1)*b-d)^2/(2*m*(m+1))):
  fi:
#case IV.4
  if  (m+1)*b<d then
    RETURN((b-a)^2/2):
  fi:
 
fi:
ERROR(`Something is wrong`):
end:
 
################################################  
 
#f122(a,b,c,d,m) finds the number of triples (i,j,i+j) s.t. a<=i<=b
#and c<=j,i+mj<=d
 
f122:=proc(a,b,c,d,m):
 
if not (0<=a and a<b and b<c and c<d) then
 ERROR(`0<=a<b<c<d`):
fi:
 
if d<=a+m*c then
 RETURN(0):
fi:
 
if a+m*c<d and d<=b+m*c then
 RETURN((d-m*c-a)^2/2):
fi:
 
if b+m*c<d then
  RETURN(((d-a-m*c)+(d-b-m*c))*(b-a)/2):
fi:
 
ERROR(`Something is wrong`):
end:
 
################################################# 
 
#f123(a,b,c,d,e,f,m) finds the number of triples (i,j,i+mj) s.t. a<=i<=b
#and c<=j<=d, and $e<=i+mj<=f
 
 
f123:=proc(a,b,c,d,e,f,m):
 
if not (0<=a and a<b and b<c and c<d and d<e and e<f) then
 ERROR(`0<=a<b<c<d<e<f`):
fi:
 
 
 
#Case I
if a+m*d<=b+m*c then
 
  #case I1 
    if e<=a+m*c then
       
      #case I1i
          if f<=a+m*c then
               RETURN(0):
          fi:
 
      #case I1ii
          if a+m*c<f and f<=a+m*d then
           RETURN((f-a-m*c)^2/2):
          fi:
 
      #case I1iii
          if a+m*d<f and f<=b+m*c then
           RETURN(((f-a-m*d)+(f-a-m*c))/2*(d-c)):
          fi:
 
      #case I1iv
          if b+m*c<f and f<=b+m*d then
           RETURN((d-c)*(b-a)-(b+m*d-f)^2/2):
          fi:
 
      #case I1v
          if b+m*d<f then
            RETURN((b-a)*(d-c)):
          fi:
 
 
  #endcase I1 
    fi:
 
  #case I2 
    if a+m*c<e and e<=a+m*d then
       
      #case I2i
          if  f<=a+m*d  then
           RETURN((f-a-m*c)^2/2-(e-a-m*c)^2/2):
          fi:
 
      #case I2ii
          if a+m*d<f and f<=b+m*c then
           RETURN((b-a)*(d-c)-(e-a-m*c)^2/2-(d-c)/2*((b-f+m*d)+(b-f+m*c))  ):
          fi:
 
      #case I2iii
          if b+m*c<f and f<=b+m*d then
           RETURN((b-a)*(d-c)-(m*d-f+b)^2/2-(e-a-m*c)^2/2  ):
          fi:
 
      #case I2iv
          if f>b+m*d then
           RETURN((b-a)*(d-c)-(e-a-m*c)^2/2  ):
          fi:
 
  #endcase I2 
    fi:
 
 
 
  #case I3 
    if a+m*d<e and e<=b+m*c then
       
      #case I3i
          if  f<=b+m*c  then
           RETURN((b-a)*(d-c)- (d-c)/2*((e-m*d-a)+(e-m*c-a))-
                                (d-c)/2*((b-f+m*c)+(b-f+m*d))     ):
          fi:
 
      #case I3ii
          if b+m*c<f and f<=b+m*d then
           RETURN((b-a)*(d-c)- (d-c)/2*((e-m*d-a)+(e-m*c-a))-
                                (b-f+m*d)^2/2     ):
          fi:
 
      #case I3iii
          if f>b+m*d then
           RETURN((b-a)*(d-c)- (d-c)/2*((e-m*d-a)+(e-m*c-a)) ):
          fi:
 
  #endcase I3 
    fi:
 
  #case I4 
    if  b+m*c<e and e<=b+m*d then
       
      #case I4i
          if  f<=b+m*d  then
            RETURN((m*d-e+b)^2/2-(m*d-f+b)^2/2):
          fi:
 
      #case I4ii
          if b+m*d<f then
            RETURN((m*d-e+b)^2/2):
          fi:
 
   fi:
  #endcase I4 
 
  #case I5 
    if  b+m*d<=e then
       RETURN(0):
    fi:
  #endcase I5
fi:
#endcase I
 
########
 
 
#Case II
if b+m*c<a+m*d then
 
  #case II1 
    if e<=a+m*c then
       
      #case II1i
          if f<=a+m*c then
               RETURN(0):
          fi:
 
      #case II1ii
          if a+m*c<f and f<=b+m*c then
           RETURN((f-a-m*c)^2/2):
          fi:
 
      #case II1iii
          if b+m*c<f and f<=a+m*d then
           RETURN(((f-b-m*c)+(f-a-m*c))/2*(b-a)):
          fi:
 
      #case II1iv
          if a+m*d<f and f<=b+m*d then
           RETURN((d-c)*(b-a)-(b+m*d-f)^2/2):
          fi:
 
      #case II1v
          if b+m*d<f then
            RETURN((b-a)*(d-c)):
          fi:
 
 
  #endcase II1 
    fi:
 
  #case II2 
    if a+m*c<e and e<=b+m*c then
       
      #case II2i
          if  f<=b+m*c  then
           RETURN((f-a-m*c)^2/2-(e-a-m*c)^2/2):
          fi:
 
      #case II2ii
          if b+m*c<f and f<=a+m*d then
           RETURN((b-a)*(d-c)-(e-a-m*c)^2/2-(b-a)/2*((m*d-f+a)+(m*d-f+b))  ):
          fi:
 
      #case II2iii
          if a+m*d<f and f<=b+m*d then
           RETURN((b-a)*(d-c)-(m*d-f+b)^2/2-(e-a-m*c)^2/2  ):
          fi:
 
      #case II2iv
          if f>b+m*d then
           RETURN((b-a)*(d-c)-(e-a-m*c)^2/2  ):
          fi:
 
  #endcase II2 
    fi:
 
 
 
  #case II3 
    if b+m*c<e and e<=a+m*d then
       
      #case II3i
          if  f<=a+m*d  then
           RETURN((b-a)*(d-c)- (b-a)/2*((e-b-m*c)+(e-a-m*c))-
                                (b-a)/2*((m*d-f+b)+(m*d-f+a))     ):
          fi:
 
      #case II3ii
          if a+m*d<f and f<=b+m*d then
           RETURN((b-a)*(d-c)- (b-a)/2*((e-b-m*c)+(e-a-m*c))-
                                (m*d-f+b)^2/2     ):
          fi:
 
      #case II3iii
          if f>b+m*d then
           RETURN((b-a)*(d-c)- (b-a)/2*((e-b-m*c)+(e-a-m*c)) ):
          fi:
 
  #endcase II3 
    fi:
 
  #case II4 
    if  a+m*d<=e and e<b+m*d then
       
      #case II4i
          if  f<b+m*d  then
            RETURN((m*d-e+b)^2/2-(m*d-f+b)^2/2):
          fi:
 
      #case II4ii
          if b+m*d<=f then
            RETURN((m*d-e+b)^2/2):
          fi:
    fi:
  #endcase II4 
 
 
  #case II5 
    if  b+m*d<=e then
       RETURN(0):
    fi:
  #endcase II5 
fi:
#endcase II
ERROR(`Something is wrong`):
end:
 
##################################################
 
Ron:=proc(resh,m)
local Efes, Ekhad,n,gu,mu,i,j,k,mu1,mu2:
 
n:=nops(resh):
Efes:=[[0,op(1,resh)]]:
Ekhad:=[]:
 
for i from 2 by 2 to nops(resh)-1 do
 Efes:=[op(Efes),[op(i,resh),op(i+1,resh)]]:
od:
 
if n mod 2=0 then
 Efes:=[op(Efes),[op(n,resh),1]]:
fi:
 
 
 
for i from 1 by 2 to nops(resh)-1 do
 Ekhad:=[op(Ekhad),[op(i,resh),op(i+1,resh)]]:
od:
 
if n mod 2=1 then
 Ekhad:=[op(Ekhad),[op(n,resh),1]]:
fi:
 
 
gu:=0:
 
for i from 1 to nops(Efes) do
mu:=op(i,Efes):
 gu:=gu+f111(op(1,mu),op(2,mu),m):
od:
 
for i from 1 to nops(Efes) do
  mu:=op(i,Efes):
 for j from i+1 to nops(Efes) do
  mu1:=op(j,Efes):
    gu:=gu+f112(op(1,mu),op(2,mu),op(1,mu1),op(2,mu1),m):
    gu:=gu+f122(op(1,mu),op(2,mu),op(1,mu1),op(2,mu1),m):
 od:
od:
 
for i from 1 to nops(Efes) do
  mu:=op(i,Efes):
 for j from i+1 to nops(Efes) do
  mu1:=op(j,Efes):
  for k from j+1 to nops(Efes) do
   mu2:=op(k,Efes):
       gu:=gu+f123(op(1,mu),op(2,mu),op(1,mu1),op(2,mu1),
                     op(1,mu2),op(2,mu2),m):
  od:
 od:
od:
 
 
 
for i from 1 to nops(Ekhad) do
mu:=op(i,Ekhad):
 gu:=gu+f111(op(1,mu),op(2,mu),m):
od:
 
for i from 1 to nops(Ekhad) do
  mu:=op(i,Ekhad):
 for j from i+1 to nops(Ekhad) do
  mu1:=op(j,Ekhad):
    gu:=gu+f112(op(1,mu),op(2,mu),op(1,mu1),op(2,mu1),m):
    gu:=gu+f122(op(1,mu),op(2,mu),op(1,mu1),op(2,mu1),m):
 od:
od:
 
for i from 1 to nops(Ekhad) do
  mu:=op(i,Ekhad):
 for j from i+1 to nops(Ekhad) do
  mu1:=op(j,Ekhad):
  for k from j+1 to nops(Ekhad) do
   mu2:=op(k,Ekhad):
       gu:=gu+f123(op(1,mu),op(2,mu),op(1,mu1),op(2,mu1),
                     op(1,mu2),op(2,mu2),m):
  od:
 od:
od:
 
gu:
end:


#m is from x+my=z
#x is the equilibrium point we want to test
conv:=proc(x,y,m)
local n,i,j,gu,k:
if nops(x)<>nops(y) then 
RETURN(`The lengths of the two lists must be equal.`):
fi:

n:=nops(x):
gu:=0:
for i from 1 to n/(m+1) do
  for j from i+1 to (n-i)/(m+1) do
  k:=i+m*j:
  gu:=gu+x[i]*(2*y[j]-x[j])+x[j]*(2*y[i]-x[j])+x[k]*(2*y[i]-x[i])
     + y[k]*(2*x[i]-y[i]) - x[k]*(2*x[j]-y[j]) - y[k]*(2*y[j]-x[j])
     + 2*(x[j]+y[j]):
  od:
od:
gu:
end: