#copyright 1998 by A. Robertson and D. Zeilberger

print(`Version of March 25, 1998`):
lprint(``):
print(`This is RON, A Maple package accompanying `):
print(`Aaron Robertrson and Doron Zeilberger's article`):
print(`"A 2-Coloring of [1,n] can have (1/22)n^2 + O(n)`):
print(`monochromatic triples, but not less!"`):
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(`100 dollars problem about finding`):
print(`katan(n):=the max of schur(a) over all 0-1 vectors of length n`):
print(`where schur(a) is the number of Schur triples: 1<=i<j<=n s.t.`):
print(`a[i]=a[j]=a[i+j]`):
print(`Graham, Rodl, and Rucinski `):
print(`proved that asymptotically, katan(n)<=n^2/18 and wondered if it can`):
print(`be improved`):
print(`In the paper by Robertson and Zeilberger (accompanying this package)`):
print(`It is proved that katan(n)=n^2/22+O(n)`):
lprint(``):
print(`contains the procedures ; B, br, checkptor, diff1, dif, Diff1, `):
print(`f2, findpol, findpol2, F, hakhi1`):
print(`hakhi2, hakhi3, katan, katand, ptor, Ron, Ronresh,tsamek`):
fi:
 
 
if nops([args])=1 and op(1,[args])=checkptor then
print(`checkptor(x,n,k) , Given a 0-1 sequence,x, of length n`):
print(` checks whether gives the signs that should be all negative`):
print(`at a local minimum`):
fi:
 
 
if nops([args])=1 and op(1,[args])=ptor then
print(`ptor(n,k) , solves the ping-pong recurrence in the`):
print(`article of Robertson and Zeilberger for n and k`):
print(`For example ptor(121,66)`):
fi:

 
 
if nops([args])=1 and op(1,[args])=diff1 then
print(`diff1(n,r): Given integers n and r`):
print(`finds diff1(F(n),r)/(2*x[r]-1) `):
fi: 
 
 
if nops([args])=1 and op(1,[args])=dif then
print(`dif(n,r): Given integers n and r`):
print(`finds the claimed expression for diff1(n,r) `):
fi: 
 
if nops([args])=1 and op(1,[args])=Diff1 then
print(`Diff1(F,r): Given an expression F of x[1], ..., x[n] on the`):
print(`Boolean cube, finds the discrete partial derivative w.r.t`):
print(`the variable x[r] `):
fi: 
 
 
if nops([args])=1 and op(1,[args])=F then
print(`F(n): Finds the quadratic form F(x[1], ..., x{n]) defined `):
print(`in the article. Here n is a positive integer`):
fi: 
 
if nops([args])=1 and op(1,[args])=Ronresh then
 
print(`Ronresh(list) gives Ron as an algebraic expression in its arguments`):
print(`in a more convenient form for later processing, i.e.`):
print(`it is given as a set of  contributions, each of which is`):
print(`a list of length 2 of the form [expr,set of linear expressions]`):
print(`where the linear expressions describe the inequalities >=0 that`):
print(`decide whether this term has to be added to Ron(resh):`):
print(`For example, try Ronresh([a]),Ronresh([a,b]`): 
fi:
 
if nops([args])=1 and op(1,[args])=findpol then
print(`findpol(f,a,a0,L,Maxdeg) finds emprically the polynomial of`):
print(` degree <=Maxdeg that fits the function f`):
print(` in the variable a near a=a0 with resolution L`):
print(`if it fails, it returns 0`):
fi: 
 
 
if nops([args])=1 and op(1,[args])=findpol2 then
print(`findpol2(f,a,b,a0,b0,L,Maxdeg) finds emprically the polynomial of`):
print(` degree <=Maxdeg that fits the function f`):
print(` in the variable a near a=a0 with resolution L`):
print(`if it fails, it returns 0`):
fi: 
 
if nops([args])=1 and op(1,[args])=katan then
print(`katan(n) finds the smallest possible number of Schur triples`):
print(`that a 0-1 vector of length n can have`):
fi:
 
 
if nops([args])=1 and op(1,[args])=bdok then
print(`bdok(resh,L) checks the validity of Ron empirically by`):
print(`comparing it with (Br(n*resh)) for L from 20 to L`):
fi:
 
if nops([args])=1 and op(1,[args])=f2 then
print(`f2(a,b) gives Ron([a,b]) in terms of the exact`):
print(`formula outputted by Ronresh, after it was cleaned up`): 
fi:
 
 
if nops([args])=1 and op(1,[args])=hakhi1 then
print(`hakhi1(L) finds the input that minimizes Ron([x]), and the`):
print(`corresponding value, in resolutions of 1/L`):
fi:
 
 
if nops([args])=1 and op(1,[args])=hakhi2 then
print(`hakhi2(L) finds the input that minimizes Ron([x,y]), x<y, and the`):
print(`corresponding value, in resolutions of 1/L`):
fi:
 
 
if nops([args])=1 and op(1,[args])=hakhi3 then
print(`hakhi3(L) finds the input that minimizes Ron([x,y,z]), x<y<z, and the`):
print(`corresponding value, in resolutions of 1/L`):
fi:
 
 
 
if nops([args])=1 and op(1,[args])=Ron then
print(`Ron(list) finds the asymptotics of schur(a)/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),type`):
print( ` Ron([4,11,10/11]) `):
fi:
 
 
if nops([args])=1 and op(1,[args])=B then
print(`B(list) finds the number of Schur triples`):
print(`of the vector [0^list[1],1^list[2],0^list[3], ...]`):
fi:
 
 
if nops([args])=1 and op(1,[args])=Br then
print(`Br(list) finds evalf(B(list)/sum(list)^2)`):
fi:
 
 
if nops([args])=1 and op(1,[args])=katanad then
print(`katanad(L) finds the list of katan(n) for n=1..L`):
fi:
 
 
if nops([args])=1 and op(1,[args])=schur then
print(`schur(a) finds the number of Schur triples in the 0-1 vector a`):
print(`in other words: the number of 1<=i<j<=n such that a[i]=a[j]=a[i+j]`):
print(`the input a is given in terms of a list`):
print(`For example schur([0,0,0,1]) should give 1`):
fi:


if nops([args])=1 and op(1,[args])=tsamek then
print(`tsamek(x) inputs a 0-1 list, x and`):
print(`finds the list [a_1,a_2, ..., ]`):
print(`such that x=0^a_1,1^a_2,0^a_3...`):
fi: 
end:
 
###################################################
################################################### 
schur:=proc(a)
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 do
 if op(i,a)=op(j,a) and op(j,a)=op(i+j,a) then
  co:=co+1:
 fi:
od:
od: 
co: 
end:


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

 
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)
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):
 
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)
local n,mu,gu:
gu:=[]:
for n from 1 to L do
print(`n=`,n):
mu:=katan(n):
print(mu):
gu:=[op(gu),mu]:
od:
gu:
end:
 
###############################################
 
B1:=proc(n)
local gu,i:
gu:=[seq(0,i=1..n)]:
schur(gu):
end:
 
###############################################
 
B:=proc(resh)
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):
end:
 
##############################################
 
Br:=proc(resh):
evalf(B(resh)/convert(resh,`+`)^2):
end:

############################################## 
 
#f111(a,b) finds the number of triples (i,j,i+j) s.t. a<=i<j<i+j<=b
f111:=proc(a,b):
 
if not (0<=a and a<b ) then
 ERROR(`0<=a<b`):
fi:
 
if b-2*a<=0 then
 RETURN(0):
else
 RETURN((b-2*a)^2/4):
fi:
end:

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

############################################# 
 
hakhi3:=proc(L)
local i,j,k,gu,mu,aluf:
 
aluf:=[]:
gu:=1:
 
for i from 1 to L-3 do
 for j from i+1 to L-2 do
  for k from j+1 to L-1 do
 mu:=Ron([i/L,j/L,k/L]):
 if mu<gu then
   gu:=mu:
   aluf:=[i/L,j/L,k/L]:
 fi:
od:
od:
od:
 
gu,aluf:
end:
 
################################################## 
 
#bdok(resh,L) checks the validity of Ron by doing Br(resh*L)
bdok:=proc(resh,L)
local gu,mu,i,k,m,n:
 
if L<20 then
ERROR(`L>=20`):
fi:
 
m:=nops(resh):
gu:=RON(resh):
 
 
 
for n from L to L do
mu:=[seq(0,i=1..trunc(op(1,resh)*n))]:
 
for k from 1 by 2 to m-2 do
mu:=[op(mu),seq(1,i=trunc(op(k,resh)*n)+1..trunc(op(k+1,resh)*n))]:
mu:=[op(mu),seq(0,i=trunc(op(k+1,resh)*n)+1..trunc(op(k+2,resh)*n))]:
od:
 
if m mod 2 =0 then
mu:=[op(mu),seq(1,i=trunc(op(m-1,resh)*n)+1..trunc(op(m,resh)*n))]:
mu:=[op(mu),seq(0,i=trunc(op(m,resh)*n)+1..n)]:
else
mu:=[op(mu),seq(1,i=trunc(op(m,resh)*n)+1..n)]:
fi:
print(evalf(schur(mu)/n^2)):
od:
 
print(`Ron(`,resh,`)=`,Ron(resh)):
end:
 
################################################### 
 
hakhim1:=proc(L,K)
local i,gu,mu,aluf,alufy,k:
 
aluf:=0:
gu:=1:
 
for i from 1 to L-1 do
 mu:=Ron([i/L]):
 if mu<gu then
   gu:=mu:
   aluf:=i/L:
 fi:
od:
 
 
for k from 2 to K do
alufy:=aluf:
for i from -(L-1) to L-1 do
 mu:=Ron([alufy+i/L^k]):
 if mu<gu then
   gu:=mu:
   aluf:=alufy+i/L^k:
 fi:
od:
od:
 
gu,aluf:
end:
 
######################################################  
 
hakhim2:=proc(L,K)
local i,j,k,gu,gu1,alufy,mu,aluf,pi,ka1,ka2,mu1,mu2:
 
aluf:=[]:
gu:=1:
 
for i from 1 to L-2 do
 for j from i+1 to L-1 do
 mu:=Ron([i/L,j/L]):
 if mu<gu then
   gu:=mu:
   aluf:=[i/L,j/L]:
 fi:
od:
od:
 
 
for k from 2 to K do
alufy:=aluf:
for i from -L to L do
for j from -L to L do
 if  0<op(1,alufy)+i/L^k and
op(1,alufy)+i/L^k<op(2,alufy)+j/L^k  and op(2,alufy)+j/L^k<1 then
 mu:=Ron([op(1,alufy)+i/L^k,op(2,alufy)+j/L^k]):
 if mu<gu then
   gu:=mu:
   aluf:=[op(1,alufy)+i/L^k,op(2,alufy)+j/L^k]:
 fi:
fi:
od:
od:
od:
 
print(`the smallest value for the Graham-Schur constant for sequences`):
print(`that change color 2 times is`):
print(gu):
print(`this occurs approximately at`):
print(aluf):
 
gu1:=convert(gu,confrac,pi):
print(`the continued fraction of the Graham-Schur number is`):
print(gu1):
print(`and its partial convergents are`):
print(pi):
 
print(`The appx. to the first component of the vector where this happens has`):
print(`a continuted fraction :`):
mu1:=convert(op(1,aluf),confrac,ka1):
print(mu1):
print(`and its partial convergents are`):
print(ka1):
 
 
print(`The appx. to the 2nd component of the vector where this happens has`):
print(`a continuted fraction :`):
mu2:=convert(op(2,aluf),confrac,ka2):
print(mu2):
print(`and its partial convergents are`):
print(ka2):
print(`To conclude, the champion is`):
print([mu1,mu2]):
gu,[mu1,mu2]:
end:
 
################################################  
 
hakhim3:=proc(L,K)
local i,j,j1,k,gu,gu1,mu,aluf,alufy,pi,ka1,ka2,ka3,  mu1,mu2,mu3 :
 
aluf:=[]:
gu:=1:
 
for i from 1 to L-3 do
 for j from i+1 to L-2 do
 for j1 from j+1 to L-1 do
 mu:=Ron([i/L,j/L,j1/L]):
 if mu<gu then
   gu:=mu:
   aluf:=[i/L,j/L,j1/L]:
 fi:
od:
od:
od:
 
for k from 2 to K do
alufy:=aluf:
for i from -L to L do
for j from -L to L do
for j1 from -L to L do
 if 0<op(1,alufy)+i/L^k and op(1,alufy)+i/L^k<op(2,alufy)+j/L^k  
     and op(2,alufy)+j/L^k<op(3,alufy)+j1/L^k  and op(3,alufy)+j1/L^k<1 then
 
       mu:=Ron([op(1,alufy)+i/L^k,op(2,alufy)+j/L^k,
               op(3,alufy)+j1/L^k ]):
 if mu<gu then
   gu:=mu:
   aluf:=[op(1,alufy)+i/L^k,op(2,alufy)+j/L^k,op(3,alufy)+j1/L^k ]:
 fi:
fi:
od:
od:
od:
od:
 
print(`the smallest value for the Graham-Schur constant for sequences`):
print(`that change color 3 times is appx.`):
print(gu):
print(`this occurs approximately at`):
print(aluf):
 
gu1:=convert(gu,confrac,pi):
print(`the continued fraction of the Graham-Schur number is`):
print(gu1):
print(`and its partial convergents are`):
print(pi):
 
print(`The appx. to the first component of the vector where this happens has`):
print(`a continuted fraction :`):
mu1:=convert(op(1,aluf),confrac,ka1):
print(mu1):
print(`and its partial convergents are`):
print(ka1):
 
print(`The appx. to the 2nd component of the vector where this happens has`):
print(`a continuted fraction :`):
mu2:=convert(op(2,aluf),confrac,ka2):
print(mu2):
print(`and its partial convergents are`):
print(ka2):
 
 
 
print(`The appx. to the 3rd component of the vector where this happens has`):
print(`a continuted fraction :`):
mu3:=convert(op(3,aluf),confrac,ka3):
print(mu3):
print(`and its partial convergents are`):
print(ka3):
 
print(`To conclude, the champion is`):
print([mu1,mu2,mu3]):
gu,[mu1,mu2,mu3]:
end:
 
###############################################  
 
hakhim4:=proc(L,K)
local i,j,j1,j2,k,gu,gu1,mu,aluf,alufy,pi,ka1,ka2,ka3,ka4,  mu1,mu2,mu3,mu4 :
 
if L<5 then
 ERROR(`L>=5`):
fi:
aluf:=[]:
gu:=1:
 
for i from 1 to L-4 do
 for j from i+1 to L-3 do
  for j1 from j+1 to L-2 do
   for j2 from j1+1 to L-1 do
 mu:=Ron([i/L,j/L,j1/L,j2/L]):
 if mu<gu then
   gu:=mu:
   aluf:=[i/L,j/L,j1/L,j2/L]:
 fi:
od:
od:
od:
od:
for k from 2 to K do
alufy:=aluf:
for i from -L to L do
for j from -L to L do
for j1 from -L to L do
for j2 from -L to L do
 if 0<op(1,alufy)+i/L^k and op(1,alufy)+i/L^k<op(2,alufy)+j/L^k  
     and op(2,alufy)+j/L^k<op(3,alufy)+j1/L^k  
       and op(3,alufy)+j1/L^k<op(4,alufy)+j2/L^k
       and op(4,alufy)+j2/L^k<1 then
 
       mu:=Ron([op(1,alufy)+i/L^k,op(2,alufy)+j/L^k,
               op(3,alufy)+j1/L^k,op(4,alufy)+j2/L^k ]):
 if mu<gu then
   gu:=mu:
   aluf:=[op(1,alufy)+i/L^k,op(2,alufy)+j/L^k,op(3,alufy)+j1/L^k,
                                             op(4,alufy)+j2/L^k ]:
 fi:
fi:
od:
od:
od:
od:
od:
print(`the smallest value for the Graham-Schur constant for sequences`):
print(`that change color 4 times is appx.`):
print(gu):
print(`this occurs approximately at`):
print(aluf):
 
gu1:=convert(gu,confrac,pi):
print(`the continued fraction of the Graham-Schur number is`):
print(gu1):
print(`and its partial convergents are`):
print(pi):
 
print(`The appx. to the first component of the vector where this happens has`):
print(`a continuted fraction :`):
mu1:=convert(op(1,aluf),confrac,ka1):
print(mu1):
print(`and its partial convergents are`):
print(ka1):
 
print(`The appx. to the 2nd component of the vector where this happens has`):
print(`a continuted fraction :`):
mu2:=convert(op(2,aluf),confrac,ka2):
print(mu2):
print(`and its partial convergents are`):
print(ka2):
 
 
 
print(`The appx. to the 3rd component of the vector where this happens has`):
print(`a continuted fraction :`):
mu3:=convert(op(3,aluf),confrac,ka3):
print(mu3):
print(`and its partial convergents are`):
print(ka3):
 
print(`The appx. to the 4th component of the vector where this happens has`):
print(`a continuted fraction :`):
mu4:=convert(op(4,aluf),confrac,ka4):
print(mu4):
print(`and its partial convergents are`):
print(ka4):
 
print(`To conclude, the champion is`):
print([mu1,mu2,mu3,mu4]):
gu,[mu1,mu2,mu3,mu4]:
end:
 
#################################################### 
  
#diff1e(resh) finds the lists of differences
 
diff1e:=proc(resh)
local gu,i:
gu:=[]:
 
for i from 1 to nops(resh)-1 do
gu:=[op(gu),op(i,resh)-op(i+1,resh)]:
od:
gu:
end:

#####################################################
 
#findDeg(f,a,a0,L,Maxdeg) finds the degree of the 
#polynomial describing f as a function of a near a0 with resolution L
#od degree<=Maxdeg. Returns -1 if not found
 
findDeg:=proc(f,a,a0,L,Maxdeg)
local i,gu,i1:
gu:=[]:
 
for i from 0 to 2*(Maxdeg+1) do
 gu:=[op(gu),f(a0+i*L)]:
od:
 
for i from 1 to Maxdeg do
gu:=diff1e(gu):
 
if gu=[seq(0,i1=1..nops(gu))] then
 RETURN(i-1):
fi:
 
od:
 
-1:
end:
 
#################################################### 
 
#findpol(f,a,a0,L,Maxdeg) finds emprically the polynomial of degree <=Maxdeg
#that fits the function f
# in the variable a near a=a0 with resolution L, 
#if it fails, it returns 0
 
findpol:=proc(f,a,a0,L,Maxdeg)
local i,deg,eq,var,b,POL:
 
deg:=findDeg(f,a,a0,L,Maxdeg+1):
 
if deg=-1 then
 RETURN(0):
fi:
 
POL:=0:
eq:={}:
var:={}:
 
for i from 0 to deg do
POL:=POL+b[i]*a^i:
var:=var union {b[i]}:
od:
 
for i from 0 to deg+1 do
eq:=eq union {subs(a=a0+i*L,POL)=f(a0+i*L)}:
od:
 
var:=solve(eq,var):
 
if var=NULL then
 ERROR(`Nothing found`):
fi:
 
POL:=subs(var,POL):
 
POL:
end:
 
################################################### 
 
#findDeg2a(f,a,b,a0,b0,L,Maxdeg) finds the degree in a of the 
#polynomial describing f as a function of a and b near (a0,b0)
# with resolution L
#of degree<=Maxdeg. Returns -1 if not found
 
findDeg2a:=proc(f,a,b,a0,b0,L,Maxdeg)
local i,gu,i1:
gu:=[]:
 
for i from 0 to 2*(Maxdeg+1) do
 gu:=[op(gu),f(a0+i*L,b0)]:
od:
 
for i from 1 to Maxdeg do
gu:=diff1e(gu):
 
if gu=[seq(0,i1=1..nops(gu))] then
 RETURN(i-1):
fi:
 
od:
 
-1:
end:
 
##################################################
 
#findDeg2b(f,a,b,a0,b0,L,Maxdeg) finds the degree in b of the 
#polynomial describing f as a function of a and b near (a0,b0)
# with resolution L
#of degree<=Maxdeg. Returns -1 if not found
 
findDeg2b:=proc(f,a,b,a0,b0,L,Maxdeg)
local i,gu,i1:
gu:=[]:
 
for i from 0 to 2*(Maxdeg+1) do
 gu:=[op(gu),f(a0,b0+i*L)]:
od:
 
for i from 1 to Maxdeg do
gu:=diff1e(gu):
 
if gu=[seq(0,i1=1..nops(gu))] then
 RETURN(i-1):
fi:
 
od:
 
-1:
end:
 
#################################################  
 
#findpol2 finds emprically the polynomial of degree <=Maxdeg
#that fits the function f
# in the variable a near a=a0 with resolution L, 
#if it fails, it returns 0
 

findpol2:=proc(f,a,b,a0,b0,L,Maxdeg)
local i1,i2,dega,degb,eq,var,c,POL,f1,f2:
 
dega:=findDeg2a(f,a,b,a0,b0,L,Maxdeg):
degb:=findDeg2b(f,a,b,a0,b0,L,Maxdeg):
 
if dega=-1 then
 RETURN(0):
fi:
 
if degb=-1 then
 RETURN(0):
fi:
 
POL:=0:
eq:={}:
var:={}:
 
for i1 from 0 to dega do
for i2 from 0 to degb do
POL:=POL+c[i1,i2]*a^i1*b^i2:
var:=var union {c[i1,i2]}:
od:
od:
 
for i1 from 0 to dega+1 do
for i2 from 0 to degb+1 do
eq:=eq union {subs({a=a0+i1*L,b=b0+i2*L},POL)=f(a0+i1*L,b0+i2*L)}:
od:
od:
var:=solve(eq,var):
 
if var=NULL then
 ERROR(`Nothing found`):
fi:
 
POL:=subs(var,POL):
 
POL:
end:

##############################################
 
graham2:=proc(a,b):
if a<=1/2 then
 1/2-b+3/4*b^2-a/4+5/4*a^2:
else
 1/2-b+3/4*b^2-a/4+1/4*a^2:
fi:
end:
 
####################################################
 
f2:=proc(a,b):
 
if not(a>=0 and b>=a and b<=1) then
   ERROR(`Wrong input`):
fi:
 
if a>=1/2 then
 1/4-b+3/4*b^2+a-a*b+1/4*a^2:
elif b<=1/2 then
1/4-b+5/4*b^2+a-2*a*b+3/4*a^2:
elif a<=1/2 and 1/2<=b and b<=1-a then
 1/4*b^2+a-2*a*b+3/4*a^2:
elif a+b>=1 and a<=1/2 then
 1/2-b+3/4*b^2-a*b+5/4*a^2:
else
 ERROR(`Too bad`):
fi:
end:
 
##################################################   
 

#Ronresh gives Ron as an algebraic expression in its arguments
#in a more convenient form for later processing, i.e.
#it is given as a set of  contribution, each of which is
#a list of length 2 of the form [expr,set of linear expressions]
#where the linear expressions describe the inequalities >=0 that
#decide whether this term has to be added to Ron(resh):
 
Ronresh1:=proc(resh)
local kakpap,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:=[]:
 
for i from 1 to nops(Efes) do
mu:=op(i,Efes):
 gu:=[op(gu),op(F111(op(1,mu),op(2,mu)))]:
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):
 kakpap:=F112(op(1,mu),op(2,mu),op(1,mu1),op(2,mu1)):
    gu:=[op(gu),op(kakpap)]:
 kakpap:=F122(op(1,mu),op(2,mu),op(1,mu1),op(2,mu1)):
    gu:=[op(gu),op(kakpap)]:
 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):
   kakpap:=F123(op(1,mu),op(2,mu),op(1,mu1),op(2,mu1),
                     op(1,mu2),op(2,mu2)):
       gu:=[op(gu),op(kakpap)]:
  od:
 od:
od:
 
 
 
for i from 1 to nops(Ekhad) do
mu:=op(i,Ekhad):
 kakpap:=F111(op(1,mu),op(2,mu)):
 gu:=[op(gu),op(kakpap)]:
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):
 kakpap:=F112(op(1,mu),op(2,mu),op(1,mu1),op(2,mu1)):
    gu:=[op(gu),op(kakpap)]:
    kakpap:= F122(op(1,mu),op(2,mu),op(1,mu1),op(2,mu1)):
    gu:=[op(gu),op(kakpap)]:
 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):
   kakpap:=F123(op(1,mu),op(2,mu),op(1,mu1),op(2,mu1),
                     op(1,mu2),op(2,mu2)):
       gu:=[op(gu),op(kakpap)]:
  od:
 od:
od:
 
gu:
end:
 
################################################# 
  
#F111(a,b) finds the number of triples (i,j,i+j) s.t. a<=i<j<i+j<=b
F111:=proc(a,b):
 
[[(b-2*a)^2/4,{b-2*a}]]:
 
end:
 
################################################### 
  
#F112(a,b,c,d) finds the number of triples (i,j,i+j) s.t. a<=i<j<=b
#and c<=i+j<=d, in terms of a,b,c,d
 
F112:=proc(a,b,c,d)
local gu:
 
 
gu:=[]:
 
  gu:=[op(gu),[(2*b-c)^2/4,{c-a-b,2*b-c,d-2*b}]]:
gu:=[op(gu),[(2*b-c)^2/4- (2*b-d)^2/4,{ c-a-b,2*b-c,2*b-d}]]:
gu:=[op(gu),[(b-a)^2/2-(c-2*a)^2/4 ,{a+b-c,c-2*a,d-2*b}]]:
gu:=[op(gu),[(b-a)^2/2-(c-2*a)^2/4-(2*b-d)^2/4 ,{a+b-c,c-2*a,d-a-b,2*b-d}]]:
gu:=[op(gu),[(d-2*a)^2/4-(c-2*a)^2/4 ,{a+b-c,c-2*a,d-2*a,a+b-d}]]:
gu:=[op(gu),[(d-2*a)^2/4 , {2*a-c,d-2*a,a+b-d}]]:
gu:=[op(gu),[(b-a)^2/2-(2*b-d)^2/4,{2*a-c,2*b-d, d-a-b}]]:
gu:=[op(gu), [(b-a)^2/2, {2*a-c,d-2*b}]]:
gu:
end:
 
############################################## 
 
#F122(a,b,c,d) finds the number of triples (i,j,i+j) s.t. a<=i<=b
#and c<=j,i+j<=d in terms of a,b,c,d
 
F122:=proc(a,b,c,d)
local gu:
 
 
gu:=[[(d-c-a)^2/2,{d-a-c,b+c-d}]]:
 
gu:=[op(gu), [((d-a-c)+(d-b-c))*(b-a)/2,{d-b-c}] ]:
 
gu:
end:
 
##################################################### 

 
#F123(a,b,c,d.e,f) finds the number of triples (i,j,i+j) s.t. a<=i<=b
#and c<=j<=d, and $e<=i+j<=f
 
 
F123:=proc(a,b,c,d,e,f)
 
local gu,
 
   guI, guI1, guI1ii, guI1iii, guI1iv, 
      guI1v, guI2, guI2i, guI2ii, guI2iii, guI2iv, guI3, guI3i, guI3ii, 
      guI3iii, guI4, guI4i, guI4ii, guII, guII1, guII1ii, guII1iii, guII1iv, 
      guII1v, guII2, guII2i, guII2ii, guII2iii, guII2iv, guII3, guII3i, 
      guII3ii, guII3iii, guII4, guII4i, guII4ii
 
:
 
gu:=[]:
 
#Case I
guI:={b+c-a-d}:
guI1:={a+c-e}:       
guI1ii:={f-a-c,a+d-f}:
gu:=[op(gu),[(f-a-c)^2/2,guI union guI1 union guI1ii]]:
guI1iii:={f-a-d,b+c-f}:
      #case I1iii
gu:=[op(gu),[((f-a-d)+(f-a-c))/2*(d-c),guI union guI1 union guI1iii]]:
 
      #case I1iv
guI1iv:={f-b-c,b+d-f}:
gu:=[op(gu),[(d-c)*(b-a)-(b+d-f)^2/2,guI union guI1 union guI1iv]]:
      #case I1v
guI1v:={f-b-d}:
 
gu:=[op(gu),[(b-a)*(d-c),guI union guI1 union guI1v]]:
 
  #endcase I1 
  #case I2 
guI2:={e-a-c,a+d-e}:
guI2i:={a+d-f}:
gu:=[op(gu),[(f-a-c)^2/2-(e-a-c)^2/2,guI union guI2 union guI2i]]:
      #case I2ii
guI2ii:={f-a-d,b+c-f}:
gu:=[op(gu),[(b-a)*(d-c)-(e-a-c)^2/2-(d-c)/2*((b-f+d)+(b-f+c))
,guI union guI2 union guI2ii]]:
      #case I2iii
guI2iii:={f-b-c,b+d-f}:
gu:=[op(gu),[(b-a)*(d-c)-(d-f+b)^2/2-(e-a-c)^2/2
,guI union guI2 union guI2iii]]:
 
      #case I2iv
guI2iv:={f-b-d}:
gu:=[op(gu),[(b-a)*(d-c)-(e-a-c)^2/2 
,guI union guI2 union guI2iv]]:
  #endcase I2 
 
  #case I3 
guI3:={e-a-d,b+c-e}:       
      #case I3i
guI3i:={b+c-f}:
 
gu:=[op(gu),[(b-a)*(d-c)- (d-c)/2*((e-d-a)+(e-c-a))-
(d-c)/2*((b-f+c)+(b-f+d)) 
,guI union guI3 union guI3i]]:
 
      #case I3ii
guI3ii:={f-b-c,b+d-f}:
 
gu:=[op(gu),[ (b-a)*(d-c)- (d-c)/2*((e-d-a)+(e-c-a))-
                                (b-f+d)^2/2
,guI union guI3 union guI3ii]]:
 
      #case I3iii
  guI3iii:={f-b-d}:
gu:=[op(gu),[ (b-a)*(d-c)- (d-c)/2*((e-d-a)+(e-c-a))
                ,guI union guI3 union guI3iii]]:
       
  #endcase I3 
 
  #case I4 
guI4:={e-b-c,b+d-e}:
 
      #case I4i
 
guI4i:={b+d-f}:
gu:=[op(gu),[ (d-e+b)^2/2-(d-f+b)^2/2
                ,guI union guI4 union guI4i]]:
 
 
      #case I4ii
 
guI4ii:={f-b-d}:
gu:=[op(gu),[ (d-e+b)^2/2
                ,guI union guI4 union guI4ii]]:
 
  #endcase I4 
 
  #case I5 
 
  #endcase I5
#endcase I
 
#Case II
 
guII:={a+d-b-c}:
  #case II1 
 
guII1:={a+c-e}:       
      #case II1i
 
      #case II1ii
 
guII1ii:={f-a-c,b+c-f}:
gu:=[op(gu),[ (f-a-c)^2/2
                ,guII union guII1 union guII1ii]]:
 
      #case II1iii
 
guII1iii:={f-b-c,a+d-f}:
gu:=[op(gu),[ ((f-b-c)+(f-a-c))/2*(b-a) 
                ,guII union guII1 union guII1iii]]:
 
 
      #case II1iv
guII1iv:={f-b-c,b+d-f}:
 
gu:=[op(gu),[ (d-c)*(b-a)-(b+d-f)^2/2
                ,guII union guII1 union guII1iv]]:
 
      #case II1v
guII1v:={f-b-d}:
 
gu:=[op(gu),[(b-a)*(d-c)
                ,guII union guII1 union guII1v]]:
 
  #endcase II1 
 
  #case II2 
guII2:={e-a-c,b+c-e}:
 
       
      #case II2i
 
guII2i:={b+c-f}:
gu:=[op(gu),[(f-a-c)^2/2-(e-a-c)^2/2
                ,guII union guII2 union guII2i]]:
 
 
      #case II2ii
 
guII2ii:={f-b-c,a+d-f}:
gu:=[op(gu),[(b-a)*(d-c)-(e-a-c)^2/2-(b-a)/2*((d-f+a)+(d-f+b))
                ,guII union guII2 union guII2ii]]:
 
 
      #case II2iii
 
 
guII2iii:={f-a-d,b+d-f}:
gu:=[op(gu),[(b-a)*(d-c)-(d-f+b)^2/2-(e-a-c)^2/2
                ,guII union guII2 union guII2iii]]:
 
 
      #case II2iv
 
guII2iv:={f-b-d}:
gu:=[op(gu),[(b-a)*(d-c)-(e-a-c)^2/2
                ,guII union guII2 union guII2iv]]:
 
 
  #endcase II2 
 
 
 
  #case II3 
 
guII3:={e-b-c,a+d-e}:
       
      #case II3i
 
guII3i:={a+d-f}:
gu:=[op(gu),[(b-a)*(d-c)- (b-a)/2*((e-b-c)+(e-a-c))-
                                (b-a)/2*((d-f+b)+(d-f+a)) 
                ,guII union guII3 union guII3i]]:
 
      #case II3ii
 
guII3ii:={f-a-d,b+d-f}:
gu:=[op(gu),[(b-a)*(d-c)- (b-a)/2*((e-b-c)+(e-a-c))-
                                (d-f+b)^2/2  
                ,guII union guII3 union guII3ii]]:
 
 
      #case II3iii
 
guII3iii:={f-b-d}:
gu:=[op(gu),[(b-a)*(d-c)- (b-a)/2*((e-b-c)+(e-a-c))
                ,guII union guII3 union guII3iii]]:
 
 
  #endcase II3 
 
 
  #case II4 
 
guII4:={e-a-d,b+d-e}:
       
      #case II4i
 
guII4i:={b+d-f}:
gu:=[op(gu),[(d-e+b)^2/2-(d-f+b)^2/2
                ,guII union guII4 union guII4i]]:
 
 
      #case II4ii
 
guII4ii:={f-b-d}:
 
gu:=[op(gu),[(d-e+b)^2/2
                ,guII union guII4 union guII4ii]]:
 
 
  #endcase II4 
 
 
  #case II5 
  #endcase II5 
 
#endcase II
 
gu:
end:
 
##############################################
 
tov:=proc(resh1)
local i:
 
 
for i from 1 to nops(resh1) do
 if op(i,resh1)<0 then
  RETURN(0):
 fi:
od:
1:
end:
 
####################################################  
 
#IsZero(kvu,resh): given a set kvu, and a list of variables resh
#[a1,...,an] decides whether or not it contains nonsenses like
#-ai,ai-aj,ai-1 1<=i<j<=n
 
IsZero:=proc(kvu,resh)
local i,j,n:
 
n:=nops(resh):
 
for i from 1 to n do
 if member(-op(i,resh),kvu) then
  RETURN(1):
fi:
 
 if member(op(i,resh)-1,kvu) then
      RETURN(1):
 fi:
od:
 
for i from 1 to n do
 for j from i+1 to n do
   if member(op(i,resh)-op(j,resh),kvu) then
    RETURN(1):
   fi:
 od:
od:
 
0:
end:
 
#################################################
 
#nake(kvu,resh): given a set kvu, and a list of variables resh
#[a1,...,an] deletes from kvu all redundant conditions
#ai,aj-aj,1-ai 1<=i<j<=n
 
nake:=proc(kvu,resh)
local i,j,n,kvu1:
 
kvu1:=kvu:
 
n:=nops(resh):
 
for i from 1 to n do
 if member(op(i,resh),kvu) then
 kvu1:=kvu1 minus {op(i,resh)}:
fi:
 
 if member(1-op(i,resh),kvu1) then
 kvu1:=kvu1 minus {1-op(i,resh)}:
 fi:
od:
 
for i from 1 to n do
 for j from i+1 to n do
   if member(op(j,resh)-op(i,resh),kvu) then
 kvu1:=kvu1 minus {op(j,resh)-op(i,resh)}:
   fi:
 od:
od:
 
kvu1:
end:
 
############################################
 
Ronresh:=proc(resh)
local gu,i,mu,gu1:
gu:=Ronresh1(resh):
 
mu:=[]:
 
for i from 1 to nops(gu) do
gu1:=op(i,gu):
 
if IsZero(op(2,gu1),resh)=0 then
 mu:=[op(mu),[op(1,gu1),nake(op(2,gu1),resh)]]:
fi:
od:
mu:
end:
 
############################################ 
 
Ronreshk:=proc(resh)
local mu,gu,gu1,gu1b,i:
gu:=Ronresh(resh):
 
mu:=[]:
 
for i from 1 to nops(gu) do
gu1:=op(i,gu): 
gu1b:=op(2,gu1):
if Fea(gu1b,resh) then
mu:=[op(mu),gu1]:
fi:
od:
mu:
end:

#################################################
 
Fea:=proc(kvu,resh)
local kvu1,kvu2,i,Tamid,n:
n:=nops(resh):
Tamid:={1-op(n,resh)}:
 
for i from 2 to n do
Tamid:=Tamid union {op(i,resh)-op(i-1,resh)}:
od:
 
kvu1:=kvu union Tamid:
 
kvu2:={}:
 
for i from 1 to nops(kvu1) do
kvu2:=kvu2 union {op(i,kvu1)>=1/10^10}:
od:
feasible(kvu2,NONNEGATIVE):
 
end:
 
################################################## 

dif:=proc(n,r)
local gu,i:
 
gu:=0:
for i from 1  to n do
gu:=gu+x[i]:
od:
 
for i from 1  to n-r do
gu:=gu+x[i]:
od:
 
gu:=gu- (n-2-trunc(r/2)+H(2*r-n-1))-x[r]:
 
if r mod 2 =0 then
gu:=gu-x[r/2]:
fi:
 
if 2*r<=n then
gu:=gu-x[r]-x[2*r]:
fi:
 
gu:
end:
 
################################################ 
 
H:=proc(a): 
if a>=0 then
 1:
else
0:
fi:
end:

############################################### 
 
F:=proc(n)
local i,j,gu:
 
gu:=0:
for i from 1 to n do
for j from i+1 to n-i do
gu:=gu+x[i]*x[j]*x[i+j]+
(1-x[i])*(1-x[j])*(1-x[i+j]):
od:
od:
gu:
end:
 
##################################################
 
diff1:=proc(n,r):
sort(factor(F(n)-subs(x[r]=1-x[r],F(n)))/(2*x[r]-1)):
end:
 
##################################################
 
Diff1:=proc(F,r):
expand(F-subs(x[r]=1-x[r],F)):
end:
 
################################################## 
 
g:=proc(n,r)
local i,gu:
 
gu:=0:
 
for i from 1 to n do
gu:=gu+x[i]:
od:
 
 
for i from 1 to n-r do
gu:=gu+x[i]:
od:
 
gu:=gu-n+r/2:
(x[r]-(1-x[r]))*gu:
end:
 
################################################
 
kapi:=proc(x)
local gu,mu,n,k,r:
gu:=[]:
n:=nops(x):
 
 
mu:=convert(x,`+`):
k:=mu:
gu:=[]:
for r from 0 to n-1 do
gu:=[op(gu),k-n+r/2+mu]:
mu:=mu-op(n-r,x):
od:
gu:
end:
 
###############################################

Fan2:=proc(x)
if x>=1 then
RETURN(0):
else
RETURN(1):
fi:
end:

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

Fan3:=proc(x)
if x>=1/2 then
RETURN(0):
else
RETURN(1):
fi:
end:
 
################################################## 
 
checkptor:=proc(x,n,k)
local r,j,lu1,lu2:
lu1:=[]:
for r from 1 to n do
lu1:=[op(lu1),sign(-1/2^20+(2*x[r]-1)*(2*k-n+r/2
-convert([op(n-r+1..n,x)],`+`)))]:
od:
 
lu2:=[]:
for r from 1 to n do
lu2:=[op(lu2),sign(-1/2^20+(2*x[r]-1)*(k-n+r/2
+convert([op(1..n-r,x)],`+`)))]:
od:
lu2,lu1:
 
end:
 
################################################# 
 
ptor:=proc(n,k)
local x,r,j:
 
for r from n to trunc((n+1)/2)+1 by -1 do
x[r]:=Fan3(k-n+trunc(r/2)+ convert([seq(x[j],j=1..n-r)],`+`) ):
x[n-r+1]:=Fan2(2*k-n-1/2+trunc((n-r+1)/2)-convert([seq(x[j],j=r..n)],`+`) ):
od:

if n mod 2=1 then
r:=trunc(n/2):
x[n-r]:=Fan3(k-n+trunc((n-r)/2)+ convert([seq(x[j],j=1..r)],`+`)):
fi:
 
[seq(x[r],r=1..n)]:
 
end:

################################################
 
kama:=proc(x,r)
local i:
for i from 1 to nops(x) while x[i]=r do
od:
i-1:
end:
 
#################################################
 

tsamek:=proc(x)
local nu,i,lu,x1,mi:
nu:=nops(x):
lu:=[]:
x1:=x:
while x1<>[] do
mi:=kama(x1,0):
lu:=[op(lu),mi]:
x1:=[op(mi+1..nops(x1),x1)]:
mi:=kama(x1,1):
lu:=[op(lu),mi]:
x1:=[op(mi+1..nops(x1),x1)]:
od:
if op(nops(lu),lu)=0 then
lu:=[op(1..nops(lu)-1,lu)]:
fi:
lu:
end:
 
#################################################