# restart; read `d:/Ake/Schur/SchurProblem.txt`; with(ListTools): print(`Welcome to SchurProblem, a Maple package written by`): print(`Thotsaporn (Aek) Thanatipanonda with help from his`): print(`lovely twin brother A. First Version: June, 27 2007 `): print(): print(`It accompanies the article: `): print(` On the monochromatic Schur Triples type problem `); print(` by Thotsaporn (Aek) Thanatipanonda`): print(`available from authors' website`): print(` http://www.math.rutgers.edu/~thot `): Help := proc() if args = NULL then print(` Part 1: LowerBound `); print(` Part 2: LowerBound2, Make`); print(` Part 31: minAllST, minAllRecursiveST `); print(` Part 32: Ord, AllTrip, EachTrip, Trip1, Trip2, Trip3, Trip4`); print(` Part 33: Zeil, RunColor, FindLength, Case1, Case2 `); print(` Case3, Case4, FindMin, FindT1T2, FindCase3 `); print(` ` ); print(` For more info: type 'Help(Part1)' or 'Help(Part2)'` ); print(` 'Help(Part31)' or 'Help(Part32)' or 'Help(Part33)' `); elif args = Part1 then print(` Function: LowerBound(k,C);`); print(` Input: the number of interval k, `); print(` list of types of upper bound C of T_{i,k-i+1} `); print(` Output: the lower bound of M(n),`); print(` the upper bound of Q, the optimal solution of Q`); print(` Try " LowerBound(2, [1]); "`); print(` Or try " LowerBound(3, [1,0]); "`); elif args = Part2 then print(` Function: LowerBound2(k,C1,C2,a);`); print(` Input: the number of interval k, `); print(` list of types of upper bound C1 and C2 of T_{i,k-i+1} `); print(` number a in equation x+ay=z`); print(` Output: the lower bound of M(n),`); print(` the upper bound of Q, the optimal solution of Q`); print(` Try " LowerBound2(4, [0,0,1,0], [1,0,0,1], 2); "`); print(` Or try " LowerBound2(4, [0,1,1,0], [1,0,0,1],3); "`); elif args = Part31 then print(` The main function of Part31: minAllSt(n,r) `); print(` Input: length of interval n, number of color r`); print(` Output: the r-coloring of all the interval of length n`); print(` that has the least number of monochromatic Schur triples`); print(` Try: 'minAllST(20,2)' `); elif args = Part32 then print(` The main function of Part32: Ord(C,L,n) `); print(` Input: the list of coloring, the list of length`); print(` corresponding to each color in C, symbol n`); print(` Output: the number of monochromatic Schur triples`); print(` Try: 'Ord([1,2,1],[1,3/2,1/4],n)' `); elif args = Part33 then print(` The main function of Part33: Zeil(r) `); print(` Input: number of color r`); print(` Output: the coloring with length of each coloring`); print(` and the number of triples of order n^2 that`); print(` obtained from Greedy Algorithm`); print(` Try: 'Zeil(2)' `); fi: end: ###################################################### # First part: the calculus part to find the # lower bound of monochromatic tripels # x+y = z ###################################################### # In C :: 1 = Area Triangle: 0 = number of R[i]*B[j] LowerBound := proc(k, C) local i, j, eq, M, R; if nops(C) <> ceil(k/2) then ERROR( nops(C) <> ceil(k/2)); elif {op(C)}minus{0,1} <> {} then ERROR(); fi: eq := add(R[i],i=1..k)*(1-add(R[i],i=1..k))+ (add(add(R[i]*(1/k-R[j])+R[j]*(1/k-R[i]) ,j=1..k-i), i = 1..k))/2 +add((1-C[i])*(R[i]*(1/k-R[k-i+1]) +R[k-i+1]*(1/k-R[i])),i=1..floor(k/2)) +add(C[i]*1^2/2/k^2,i=1..floor(k/2)); if ceil(k/2) <> floor(k/2) then i := ceil(k/2); eq := eq + (1-C[i])*(R[i]*(1/k-R[k-i+1]) +R[k-i+1]*(1/k-R[i]))/2+C[i]*1^2/2/k^2/2; fi: #print(eq); M := maximize(eq,seq(op({0<=R[i], R[i]<=1/k}), i =1..k),location); return( [(1/2-M[1])/2, M[1], M[2]]); end: ################################################################## # Second part: the calculus part to find # the lower bound of monochromatic # triples x+ay = z # ################################################################## #C1 and C2 are list of 0 or 1 #0 = Red[i]*Blue[j], 1 = area of intersection LowerBound2 := proc(k, C1, C2, a) option remember; local i, j, res, eq, M, R,x; global max1,pos, L; if nops(C1) <> k then ERROR( nops(C1) <> k); elif nops(C2) <> k then ERROR( nops(C2) <> k); elif {op(C1)} union {op(C2)} minus{0,1} <> {} then ERROR(); fi: res := k mod a; eq := add(R[i],i=1..k)*(1-add(R[i],i=1..k))/a +add(add(Make(i,j,k,R),j=1..k-a*i),i=1..floor(k/a)) +add(add(Make(i,j,k,R),j=a*i+1..k),i=1..floor(k/a)) +add((1-C1[i])*(R[ceil(i/a)]*(1/k-R[i]) +R[i]*(1/k-R[ceil(i/a)])),i=1..k) +add((1-C2[i])*(R[ceil(i/a)]*(1/k-R[k-i+1]) +R[k-i+1]*(1/k-R[ceil(i/a)])),i=1..k) +add(add(C1[a*(i-1)+j]*(1/2)*(1/k)*(j-1+j)/(a*k),j=1..a),i=1..floor(k/a)) +add(add(C2[a*(i-1)+j]*(1/2)*(1/k)*(j-1+j)/(a*k),j=1..a),i=1..floor(k/a)) +add(C1[a*(ceil(k/a)-1)+j]*(1/2)*(1/k)*(j-1+j)/(a*k),j=1..res) +add(C2[a*(ceil(k/a)-1)+j]*(1/2)*(1/k)*(j-1+j)/(a*k),j=1..res); #return(eq); eq := expand(subs({seq(R[i]=(1+x[i])/(2*k),i=1..k)},eq)); # shift and rescale R[i] #return(eq); M := maximize(eq,seq(op({-1<=x[i], x[i]<=1}), i =1..k),location); return( [(1/a-M[1])/2, M[1], M[2]]); # the first one give lower bound # M[1] is an upper bound of Q # M[2] is the set of solution that give upper bound end: Make := proc(a,b,k,R); return( R[a]*(1/k-R[b])+ R[b]*(1/k-R[a])); end: ################################################################## # Third part: Greedy Algorithm # ################################################################## # Part 31) # Main Function - Min All Schur Triple # Input: length of interval, number of color # Output: minimum value, set of minimum interval # minAllST := proc(n, r) option remember; global minValue, minInterval; minValue := n*(n-1)/2; minInterval := []; minAllRecursiveST(n, r, [1], 0); return [minValue, minInterval]; end: # # Min All Recursive Schur Triple # minAllRecursiveST := proc(n, r, interval, count) global minValue, minInterval; local i, size, count2; size := nops(interval); count2 := count; for i from 1 to size/2 do if interval[i] = interval[size] and interval[i] = interval[size-i] then count2 := count2+1; fi; od; if count2 > minValue then return; elif size = n then if count2 < minValue then minValue := count2; minInterval := [interval]; elif count2 = minValue then minInterval := [op(minInterval), interval]; fi; else for i from 1 to r do minAllRecursiveST(n, r, [op(interval), i], count2); od; fi; end: ############################################################################# ############################################################################# ## Part 32) ## Count Triple and give order. # Ord # input : list of color and length of each color # output : the order of number of Schur Triples. # Try : Ord([1,2,1],[4,6,1],n); Ord := proc(color,length,n) local nT, L; nT := AllTrip(color, length); L := add(length[i] , i = 1..nops(length)); return simplify(subs( k = n/L, nT)); end: # input: color := [1,2,1] , length := [4,6,1] AllTrip := proc(color, length) local i,j,k,numcol,Atrip, temp; numcol := max(op(color)); Atrip := 0; for i from 1 to numcol do temp := []; for j from 1 to nops(color) do if color[j] = i then temp := [op(temp), [add(length[k], k = 1..j-1),add(length[k], k = 1..j)]]; fi: od: Atrip := Atrip + EachTrip(temp); od: return(Atrip); end: # EachTrip :: # input list of interval # output number of triple EachTrip := proc(AA) local i,j,m,Etrip; Etrip := 0; for i from 1 to nops(AA) do Etrip := Etrip + Trip1(AA[i]); for j from i+1 to nops(AA) do Etrip := Etrip + Trip2(AA[i],AA[j])+ Trip3(AA[i],AA[j]); for m from j+1 to nops(AA) do Etrip := Etrip + Trip4(AA[i],AA[j],AA[m]); od: od: od: return simplify(Etrip); end: # Trip1: three points in T1 Trip1 := proc(T1) local a1,a2; a1 := T1[1]; a2 := T1[2]; return (max(a2-2*a1,0)^2/4*k^2); end: # Trip2 : one point in T1 , two point in T2. Trip2 := proc(T1,T2) local a1,a2,b1,b2,L; a1 := T1[1]; a2 := T1[2]; b1 := T2[1]; b2 := T2[2]; if not(a1>=0) or not(a1 <= a2) or not(a2 <= b1) or not(b1 <= b2) then ERROR(BadInput); fi: L := a2-a1; if b2 > b1 + a2 then return ((b2-a2-b1)*L + L^2/2)*k^2; else return max(b2-a1-b1,0)^2/2*k^2; fi: end: # Trip3 : two points in T1 , one point in T2. Trip3 := proc(T1,T2) local a1,a2,b1,b2,x,y,z,w; a1 := T1[1]; a2 := T1[2]; b1 := T2[1]; b2 := T2[2]; if not(a1>=0) or not(a1 <= a2) or not(a2 <= b1) or not(b1 <= b2) then ERROR(BadInput); fi: x := b1-a2; y := b1/2; z := b2-a2; w := b2/2; if y <= x or 2*a1 >= b2 then return(0); elif z < y then if a1 < z then return (((z-max(a1,x))*(max(a1-x,0)+z-x)/2 + (y-z)*(z-x)+ (z-x)^2/4)*k^2); elif a1 < y then return( ((y-a1)*(z-x)+ (z-x)^2/4)*k^2); elif a1 < w then return( (1/2*(w-a1)*2*(w-a1))*k^2); fi: elif z < w then if a1 < y then return ( ((y-max(a1,x))*(max(a1-x,0)+y-x)/2 +1/2*(z-y)*(3*y-2*x-z) + (2*y-x-z)^2/4)*k^2); elif a1 < z then return ( (1/2*(z-a1)*(2*y-x-z+2*y-x-a1)+(2*y-x-z)^2/4)*k^2); elif a1 < w then return ( ((w-a1)*2*(w-a1)/2)*k^2); fi: elif z >= w then if a1 < y then return (((y-max(a1,x))*(max(a1-x,0)+y-x)/2+(y-x)^2/2)*k^2); elif a1 < w then return ((2*y-x-a1)^2/2*k^2); fi: else ERROR("No case in Trip3."); fi: end: # Trip4 : one point in T1, one point in T2, one point in T3. Trip4 := proc(T1,T2,T3) local a1,a2,b1,b2,c1,c2,x,y,z,w; a1 := T1[1]; a2 := T1[2]; b1 := T2[1]; b2 := T2[2]; c1 := T3[1]; c2 := T3[2]; if not(a1>=0) or not(a1 <= a2) or not(a2 <= b1) or not(b1 <= b2) or not(b2 <= c1) or not(c1 <= c2) then ERROR(BadInput); fi: x := c1-a2; y := min(c1-a1,c2-a2); z := max(c1-a1,c2-a2); w := c2-a1; if w <= b1 or x >= b2 then return(0); elif b1 < y then if b2 < y then return(((b2-max(b1,x))*(max(b1-x,0)+b2-x)/2)*k^2); elif b2 < z then return(((y-max(b1,x))*(max(b1-x,0)+y-x)/2 + (y-x)*(b2-y))*k^2); else return(((y-max(b1,x))*(max(b1-x,0)+y-x)/2 + (y-x)*(z-y) + 1/2*(min(b2,w)-z)*(y-x+max(w-b2,0)))*k^2); fi: elif b1 < z then if b2 < z then return( (b2-b1)*(y-x)*k^2); else return( ((z-b1)*(y-x) + 1/2*(min(b2,w)-z)*(y-x+max(w-b2,0)))*k^2); fi: elif b1 < w then return( (1/2*(min(b2,w)-b1)*(w-b1+max(w-b2,0)))*k^2); else ERROR("No case in Trip4"); fi: end: ############################################################################ ############################################################################ # Part 33) # Use calculus and greedy algorithm to find the coloring that give the # minimum number of Schur Triples. # Zeil # input : color # output : the minimum coloring, length of each coloring, the order. Zeil := proc(r) local color , length , temp; color := [1]; length := [1]; temp := RunColor(color,length,r); while temp <> [color,length] do color := temp[1]; length := temp[2]; temp := RunColor(color,length,r); od: return(color, length, Ord(color, length,n)); end: # RunColor # input : color, length, r # output : the new minimum coloring, new length of each coloring. # of the smallest possible color (according to greedy algorithm). RunColor := proc(color, length, r) local i, temp; for i from 1 to r do temp := FindLength(color, length, i); if temp <> 0 then return([[op(color),i], [op(length),temp]]); fi: od: return([color,length]); end: # FindLength # input : color, length, index # output : the new optimal length (attaching at the end) # of the index coloring. FindLength := proc(color, length, index) local n,ret, more, newL; if nops(color) = 1 and index = 1 then return 0; fi: more := Case1(color, length, index, var) + Case2(color, length, index,var) + Case3(color, length, index,var) + Case4(color, length, index,var); newL := FindMin(color, length, index,more , var); ############################################## #check output, I check the output in FindMin;# ############################################## return(newL); end: ############################################################################## ############################################################################## # Case 1 Case1 := proc(color, length, index, X) local long; long := add(length[i], i = 1..nops(length)); return( piecewise(X > long, (X-long)^2/4*k^2)); end: #Case 2 Case2 := proc(color, length, index, X) local a1,a2,b1,L; if not(member(index, {op(color)})) then return 0; fi: a1 := add(length[i], i = 1..Search(index, color)-1); a2 := add(length[i], i = 1..Search(index, color)); b1 := add(length[i], i = 1..nops(length)); L := a2-a1; return(piecewise( X > a2,((X-a2)*L + L^2/2)*k^2, X >= 0, max(X-a1,0)^2/2*k^2)); end: #Case 3 Case3 := proc(color, length, index, var) local a1,a2,b1,b2 ,x,y,z,w, endPos, listTT; listTT := FindCase3(color, length, index); if listTT = [] then return(0); fi: endPos := add(length[i], i = 1..nops(length)); a1 := listTT[1]; a2 := listTT[2]; b1 := endPos; b2 := endPos+var; x := b1-a2; y := b1/2; z := b2-a2; w := b2/2; if y <= x then return(0); fi: return piecewise( a1 <= x and z <= y, (1/2*(z-x)^2+(y-z)*(z-x)+ (z-x)^2/4)*k^2, a1 <= z and z <= y, (1/2*(z-a1)*(z-x+a1-x)+(y-z)*(z-x)+ (z-x)^2/4)*k^2, a1 <= y and z <= y, ((y-a1)*(z-x)+(z-x)^2/4)*k^2, a1 <= w and z <= y, ((w-a1)*2*(w-a1)/2)*k^2, a1 <= x and z <= w, ((y-x)^2/2 + 1/2*(z-y)*(y-x+2*y-x-z)+ 1/4*(2*y-x-z)^2)*k^2, a1 <= y and z <= w, (1/2*(y-a1)*(y-x+a1-x) +1/2*(z-y)*(y-x+2*y-x-z)+1/4*(2*y-x-z)^2)*k^2, a1 <= z and z <= w, (1/2*(z-a1)*(y-x-(a1-y)+2*y-x-z)+1/4*(2*y-x-z)^2)*k^2, a1 <= w and z <= w, ((w-a1)*2*(w-a1)/2)*k^2, a1 <= x and w < z, ((y-x)^2)*k^2, a1 <= y and w < z, (1/2*(y-a1)*(y-x+a1-x)+1/2*(y-x)^2)*k^2, a1 <= w and w < z, ((y+y-x-a1)^2/2)*k^2, 0); end: #Case4 Case4 := proc(color, length, index, var) local a1,a2,b1,b2,c1,c2,x,y,z,w ,i, listTT, f,g; listTT := FindT1T2( color, length, index); f := 0 ; for i from 1 to nops(listTT) do a1 := listTT[i][1][1]; a2 := listTT[i][1][2]; b1 := listTT[i][2][1]; b2 := listTT[i][2][2]; c1 := add(length[j], j = 1..nops(length)); c2 := c1+var; x := c1-a2; y := c2-a2; z := c1-a1; w := c2-a1; g := piecewise( b2 < x , 0, b2 < y and y < z, ((b2-x)^2/2)*k^2, b2 < z and y < z, ((y-x)^2/2 + (b2-y)*(y-x))*k^2 , b2 < w and y < z, ((y-x)^2/2 + (z-y)*(y-x) + (b2-z)*((y-x)+w-b2)/2)*k^2, y < z, ((y-x)^2/2 + (z-y)*(y-x) + (w-z)^2/2)*k^2, b2 < z , ((b2-x)^2/2)*k^2, b2 < y , ((z-x)^2/2 + (b2-z)*(z-x))*k^2 , b2 < w , ((z-x)^2/2 + (y-z)*(z-x) + (b2-y)*((z-x)+w-b2)/2)*k^2 , ((z-x)^2/2 + (y-z)*(z-x) + (w-y)^2/2)*k^2 ) - piecewise( b1 < x , 0, b1 < y and y < z, ((b1-x)^2/2)*k^2, b1 < z and y < z, ((y-x)^2/2 + (b1-y)*(y-x))*k^2 , b1 < w and y < z, ((y-x)^2/2 + (z-y)*(y-x) + (b1-z)*((y-x)+w-b1)/2)*k^2, y < z, ((y-x)^2/2 + (z-y)*(y-x) + (w-z)^2/2 )*k^2, b1 < z , ((b1-x)^2/2)*k^2, b1 < y , ((z-x)^2/2 + (b1-z)*(z-x))*k^2 , b1 < w , ((z-x)^2/2 + (y-z)*(z-x) + (b1-y)*((z-x)+w-b1)/2)*k^2 , ((z-x)^2/2 + (y-z)*(z-x) + (w-y)^2/2)*k^2); f := f + g; od: return f; end: # FindT1T2 for case 4. FindT1T2 := proc(color,length,index) local i, j, listCol, endPos, ret; endPos := add(length[i], i = 1..nops(length)); listCol := []; ret := []; #make list color for i from 1 to nops(color) do if color[i] = index then listCol := [op(listCol), [add(length[m], m = 1..i-1),add(length[m], m = 1..i)]]; fi: od: for i from 1 to nops(listCol) do for j from i +1 to nops(listCol) do if ( listCol[i][1] + listCol[j][1] <= endPos and listCol[i][2] + listCol[j][2] >= endPos) then ret := [op(ret),[listCol[i],listCol[j]]]; fi: od: od: return ret; end: FindCase3 := proc(color, length, index) local i, endPos ; endPos := add(length[i], i = 1..nops(length)); i := 1; while i < nops(color) and add(length[j],j=1..i) < endPos/2 do i := i +1; od: if color[i] = index then return( [add(length[j],j=1..i-1) ,add(length[j],j=1..i)] ); else return( [] ); fi: end: #FindMin FindMin := proc(color, length, index,more, var) local n, long, nT ,EQ, temp ; long := add(length[i], i = 1..nops(length)); nT := AllTrip(color,length) + more; assume(n >= 0); EQ := simplify(subs(k = n/(long + var),nT)) ; temp := minimize( EQ , var >= 0, location); #check output if (evalb( temp[1] = Ord( [op(color),index],[op(length),rhs(temp[2][1][1][1])],n))) = false then print("Counting triples meet something not expected"); ERROR(); fi: return rhs(temp[2][1][1][1]); end: