#
# Maple spreadsheet used to find coefficient for
# a block pattern coloring for monochromatic
# constellations. A supplement to the paper
# "Finding patterns avoiding many monochromatic
# constellations" by Steve Butler, Kevin Costello
# and Ron Graham
#
# Written by Steve Butler
# last updated on 10 June 2009
#
# load the following procedures
# (examples on how to use this are given below)
#
# LinearAlgebra package is only used to solve
# linear system of equations, everything else
# is built into the below procedures
with(LinearAlgebra):
#
# computes coefficient given constellation and
# subdivision pattern
Constellation:=proc(Cpat,Spat) # compute the coefficient
global CCoeff,preCCoeff,perturb,M,Q,Spatnew,blocks;
local i,j;
CCoeff:=0;
if (perturb=true) then # if perturbing we need to load some matrices
M:=Matrix(nops(Spat),nops(Spat));
Q:=Vector(nops(Spat));
Q[nops(Spat)]:=1;
end if;
for i from 1 by 2 while i<nops(Spat) do # start the recursions for red
for j from 1 by 2 while j<nops(Spat) do
RECURSEP([[Spat[i],Spat[j],1,j],[Spat[i+1],Spat[j],nops(Cpat),i+1],[Spat[i+1],Spat[j+1],1,j+1],[Spat[i],Spat[j+1],nops(Cpat),i]],Cpat,Spat,2,0);
end do;
end do;
for i from 2 by 2 while i<nops(Spat) do
for j from 2 by 2 while j<nops(Spat) do # start the recursions for blue
RECURSEP([[Spat[i],Spat[j],1,j],[Spat[i+1],Spat[j],nops(Cpat),i+1],[Spat[i+1],Spat[j+1],1,j+1],[Spat[i],Spat[j+1],nops(Cpat),i]],Cpat,Spat,2,1);
end do;
end do;
i:=true; # test for symmetry
for j from 1 to nops(Cpat) do
if (Cpat[j]+Cpat[nops(Cpat)+1-j]<>1) then i:=false; end if;
end do;
if (i=true) then CCoeff:=CCoeff/2; end if;
i:=1; # next get common denominator
for j from 2 to (nops(Cpat)-1) do
i:=lcm(i,denom(Cpat[j]));
end do;
if (perturb=true) then
preCCoeff:=CCoeff/i;
printf("The pattern before perturbing gives coefficient:");
print(preCCoeff);
perturb:=false;
Spatnew:=Spat;
for i from 1 to nops(Spat) do M[1,i]:=0; M[nops(Spat),i]:=0; end do; # finish editing M
M[1,1]:=1;
M[nops(Spat),nops(Spat)]:=1;
Spatnew:=convert(LinearSolve(M,Q),list); # now solve the linear system of equations
i:=true;
for j from 1 to (nops(Spatnew)-1) do # check to see if answer is reasonable
if (Spatnew[j]>Spatnew[j+1]) then i:=false; end if;
end do;
if (i=false) then printf("The pattern found by perturbing is inconsistent.");
else
i:=1; blocks:=Vector(nops(Spatnew)-1);
for j from 1 to (nops(Spatnew)-1) do # find new block structure
blocks[j]:=Spatnew[j+1]-Spatnew[j];
i:=lcm(i,denom(blocks[j]));
end do;
for j from 1 to (nops(Spatnew)-1) do
blocks[j]:=i*blocks[j];
end do;
blocks:=convert(blocks,list);
printf("Perturbing suggests the following pattern:");
print(blocks);
printf("The pattern after perturbing gives coefficient:");
Constellation(Cpat,Spatnew);
perturb:=true;
printf("This gives a difference (original)-(perturbed)");
preCCoeff-CCoeff;
end if;
else
CCoeff:=CCoeff/i;
print(CCoeff);
end if;
end proc:
#
# find the area of a triangle with points A, B, C
AREAT:=proc(A,B,C)
local AA, BB, CC;
if (A[1]<=B[1] and B[1]<=C[1]) then
AA:=A; BB:=B; CC:=C;
elif (A[1]<=C[1] and C[1]<=B[1]) then
AA:=A; BB:=C; CC:=B;
elif (B[1]<=A[1] and A[1]<=C[1]) then
AA:=B; BB:=A; CC:=C;
elif (B[1]<=C[1] and C[1]<=A[1]) then
AA:=B; BB:=C; CC:=A;
elif (C[1]<=A[1] and A[1]<=B[1]) then
AA:=C; BB:=A; CC:=B;
else
AA:=C; BB:=B; CC:=A;
end if;
if (AA[2]<=BB[2] and BB[2]<=CC[2]) then
if ((BB[2]-AA[2])*(CC[1]-AA[1])>(CC[2]-AA[2])*(BB[1]-AA[1])) then
(CC[1]-AA[1])*(CC[2]-AA[2])/2-(BB[1]-AA[1])*(CC[2]-BB[2])-(CC[1]-BB[1])*(CC[2]-BB[2])/2-(BB[1]-AA[1])*(BB[2]-AA[2])/2;
else
(CC[1]-AA[1])*(CC[2]-AA[2])/2-(CC[1]-BB[1])*(BB[2]-AA[2])-(CC[1]-BB[1])*(CC[2]-BB[2])/2-(BB[1]-AA[1])*(BB[2]-AA[2])/2;
end if;
elif (AA[2]<=CC[2] and CC[2]<=BB[2]) then
(CC[1]-AA[1])*(BB[2]-AA[2])-(BB[1]-AA[1])*(BB[2]-AA[2])/2-(CC[1]-BB[1])*(BB[2]-CC[2])/2-(CC[1]-AA[1])*(CC[2]-AA[2])/2;
elif (BB[2]<=AA[2] and AA[2]<=CC[2]) then
(CC[1]-AA[1])*(CC[2]-BB[2])-(BB[1]-AA[1])*(AA[2]-BB[2])/2-(CC[1]-BB[1])*(CC[2]-BB[2])/2-(CC[1]-AA[1])*(CC[2]-AA[2])/2;
elif (BB[2]<=CC[2] and CC[2]<=AA[2]) then
(CC[1]-AA[1])*(AA[2]-BB[2])-(BB[1]-AA[1])*(AA[2]-BB[2])/2-(CC[1]-BB[1])*(CC[2]-BB[2])/2-(CC[1]-AA[1])*(AA[2]-CC[2])/2;
elif (CC[2]<=AA[2] and AA[2]<=BB[2]) then
(CC[1]-AA[1])*(BB[2]-CC[2])-(BB[1]-AA[1])*(BB[2]-AA[2])/2-(CC[1]-BB[1])*(BB[2]-CC[2])/2-(CC[1]-AA[1])*(AA[2]-CC[2])/2;
else
if ((AA[2]-BB[2])*(AA[1]-CC[1])>(AA[2]-CC[2])*(AA[1]-BB[1])) then
(CC[1]-AA[1])*(AA[2]-CC[2])/2-(CC[1]-BB[1])*(AA[2]-BB[2])-(CC[1]-BB[1])*(BB[2]-CC[2])/2-(BB[1]-AA[1])*(AA[2]-BB[2])/2;
else
(CC[1]-AA[1])*(AA[2]-CC[2])/2-(BB[1]-AA[1])*(BB[2]-CC[2])-(CC[1]-BB[1])*(BB[2]-CC[2])/2-(BB[1]-AA[1])*(AA[2]-BB[2])/2;
end if;
end if;
end proc:
#
# the area of (convex) polygon P
AREAP:=proc(P)
local i,n,area;
n:=nops(P);
area:=0;
for i from 2 to n-1 do
area:=area+AREAT(P[1],P[i],P[i+1]);
end do;
area;
end proc:
#
# find the intersection of line between A and B and line y=mx+b
INTERSECT:=proc(A,B,m,b)
local mm,bb,U;
if (A[1]=B[1]) then
U[1]:=A[1];
U[2]:=m*A[1]+b;
else
mm:=(B[2]-A[2])/(B[1]-A[1]);
bb:=A[2]-mm*A[1];
U[1]:=(bb-b)/(m-mm);
U[2]:=m*U[1]+b;
end if;
convert(U,list);
end proc:
#
# The heart of the worksheet. The idea is to take P
# a polygon and then we slice it off according to the
# current set of slopes that we are dealing with
# note that each point is a quadruple, the first
# two entries is the location of the point and the second
# two give slope and intercept (indirectly) of line
# on incident edge.
RECURSEP:=proc(P,Cpat,Spat,iterate,MOD)
global CCoeff;
local i,j,k,Q,R,S,INTER,numABOVE,numBELOW,DONE,ptA,ptAloc,ptAfound,ptB,ptBloc,cntr;
if (iterate=nops(Cpat)) then
CCoeff:=CCoeff+AREAP(P);
if (perturb=true) then loadM(P,Cpat,Spat,MOD); end if; # load the matrix if perturbing
else
Q:=P;
DONE:=false;
for i from 2 by 1 while (i<=nops(Spat) and DONE=false) do # see whether points of P are above or below in current line
numABOVE:=0; numBELOW:=0;
for j from 1 to nops(Q) do
if (Q[j][2]<(Spat[i]-Cpat[iterate]*Q[j][1])/(1-Cpat[iterate])) then
numBELOW:=numBELOW+1;
INTER[j]:=-1;
else
numABOVE:=numABOVE+1;
INTER[j]:=1;
end if;
end do;
if (numBELOW>0) then # we can do something
if (numABOVE=0) then #everything is between previous and current line
DONE:=true;
if ((i-MOD) mod 2=0) then
RECURSEP(Q,Cpat,Spat,iterate+1,MOD);
end if;
else # have to slice part of Q off
ptAfound:=false; # we now find how it slices in half
for k from 1 to (nops(Q)-1) do
if (INTER[k]*INTER[k+1]=-1) then
if (ptAfound=false) then
ptAfound:=true;
ptAloc:=k;
ptA:=INTERSECT(Q[k],Q[k+1],Cpat[iterate]/(Cpat[iterate]-1),Spat[i]/(1-Cpat[iterate]));
else
ptBloc:=k;
ptB:=INTERSECT(Q[k],Q[k+1],Cpat[iterate]/(Cpat[iterate]-1),Spat[i]/(1-Cpat[iterate]));
end if;
end if;
end do;
if (INTER[nops(Q)]*INTER[1]=-1) then
ptBloc:=nops(Q);
ptB:=INTERSECT(Q[1],Q[nops(Q)],Cpat[iterate]/(Cpat[iterate]-1),Spat[i]/(1-Cpat[iterate]));
end if;
unassign('R');
unassign('S');
cntr:=1;
for k from 1 to ptAloc do
R[cntr]:=Q[k];
cntr:=cntr+1;
end do;
R[cntr][1]:=ptA[1]; R[cntr][2]:=ptA[2]; R[cntr][3]:=iterate; R[cntr][4]:=i;
R[cntr]:=convert(R[cntr],list);
cntr:=cntr+1;
R[cntr][1]:=ptB[1]; R[cntr][2]:=ptB[2]; R[cntr][3]:=Q[ptBloc][3]; R[cntr][4]:=Q[ptBloc][4];
R[cntr]:=convert(R[cntr],list);
cntr:=cntr+1;
for k from (ptBloc+1) to nops(Q) do
R[cntr]:=Q[k];
cntr:=cntr+1;
end do;
cntr:=1;
for k from (ptAloc+1) to ptBloc do
S[cntr]:=Q[k];
cntr:=cntr+1;
end do;
S[cntr][1]:=ptB[1]; S[cntr][2]:=ptB[2]; S[cntr][3]:=iterate; S[cntr][4]:=i;
S[cntr]:=convert(S[cntr],list);
cntr:=cntr+1;
S[cntr][1]:=ptA[1]; S[cntr][2]:=ptA[2]; S[cntr][3]:=Q[ptAloc][3]; S[cntr][4]:=Q[ptAloc][4];
S[cntr]:=convert(S[cntr],list);
if (INTER[ptAloc+1]=1) then # S is above line keep it, recurse on R
Q:=convert(S,list);
if ((i-MOD) mod 2=0) then
RECURSEP(convert(R,list),Cpat,Spat,iterate+1,MOD);
end if;
else # S is below recurse it, keep R
Q:=convert(R,list);
if ((i-MOD) mod 2=0) then
RECURSEP(convert(S,list),Cpat,Spat,iterate+1,MOD);
end if;
end if;
end if;
end if;
end do;
end if;
end proc:
#
# the following loads the entries of M according to the polygon Pin
loadM:=proc(Pin,Cpat,Spat,MOD)
global M;
local i,n,P;
n:=nops(Pin);
for i from 1 to n do
P[i]:=Pin[i];
end do;
P[n+1]:=P[1];
P[n+2]:=P[2];
P:=convert(P,list);
for i from 2 to (nops(P)-1) do
if (P[i][3]=nops(Cpat)) then # a vertical line
if ((P[i][2]<P[i+1][2] and MOD=0) or (P[i][2]>P[i+1][2] and MOD=1)) then
M[P[i][4],P[i-1][4]]:=M[P[i][4],P[i-1][4]]-Cpat[P[i][3]]/(Cpat[P[i-1][3]]-Cpat[P[i][3]]);
M[P[i][4],P[i][4]]:=M[P[i][4],P[i][4]]+Cpat[P[i-1][3]]/(Cpat[P[i-1][3]]-Cpat[P[i][3]])+Cpat[P[i+1][3]]/(Cpat[P[i][3]]-Cpat[P[i+1][3]]);
M[P[i][4],P[i+1][4]]:=M[P[i][4],P[i+1][4]]-Cpat[P[i][3]]/(Cpat[P[i][3]]-Cpat[P[i+1][3]]);
else
M[P[i][4],P[i-1][4]]:=M[P[i][4],P[i-1][4]]+Cpat[P[i][3]]/(Cpat[P[i-1][3]]-Cpat[P[i][3]]);
M[P[i][4],P[i][4]]:=M[P[i][4],P[i][4]]-Cpat[P[i-1][3]]/(Cpat[P[i-1][3]]-Cpat[P[i][3]])-Cpat[P[i+1][3]]/(Cpat[P[i][3]]-Cpat[P[i+1][3]]);
M[P[i][4],P[i+1][4]]:=M[P[i][4],P[i+1][4]]+Cpat[P[i][3]]/(Cpat[P[i][3]]-Cpat[P[i+1][3]]);
end if;
else
if ((P[i][1]>P[i+1][1] and MOD=0) or (P[i][1]<P[i+1][1] and MOD=1)) then
M[P[i][4],P[i-1][4]]:=M[P[i][4],P[i-1][4]]+1/(Cpat[P[i][3]]-Cpat[P[i-1][3]]);
M[P[i][4],P[i][4]]:=M[P[i][4],P[i][4]]+((1-Cpat[P[i+1][3]])/(Cpat[P[i][3]]-Cpat[P[i+1][3]])+(1-Cpat[P[i-1][3]])/(Cpat[P[i-1][3]]-Cpat[P[i][3]]))/(1-Cpat[P[i][3]]);
M[P[i][4],P[i+1][4]]:=M[P[i][4],P[i+1][4]]+1/(Cpat[P[i+1][3]]-Cpat[P[i][3]]);
else
M[P[i][4],P[i-1][4]]:=M[P[i][4],P[i-1][4]]-1/(Cpat[P[i][3]]-Cpat[P[i-1][3]]);
M[P[i][4],P[i][4]]:=M[P[i][4],P[i][4]]-((1-Cpat[P[i+1][3]])/(Cpat[P[i][3]]-Cpat[P[i+1][3]])+(1-Cpat[P[i-1][3]])/(Cpat[P[i-1][3]]-Cpat[P[i][3]]))/(1-Cpat[P[i][3]]);
M[P[i][4],P[i+1][4]]:=M[P[i][4],P[i+1][4]]-1/(Cpat[P[i+1][3]]-Cpat[P[i][3]]);
end if;
end if;
end do;
end proc:
#
# The following is a simplified version of using
# Constellation, A is any realization of the constellation
# and B is the block structure.
blockConstellation:=proc(A,B)
local i,AA,BB,nA,nB,total;
nA:=nops(A);
nB:=nops(B);
for i from 1 to nA do
AA[i]:=(A[i]-A[1])/(A[nA]-A[1]);
end do;
total:=0;
for i from 1 to nB do
total:=total+B[i];
end do;
BB[1]:=0;
for i from 1 to nB do
BB[i+1]:=BB[i]+B[i]/total;
end do;
Constellation(convert(AA,list),convert(BB,list));
end proc:#
# the main function is blockConstellation. There are two
# inputs, one a constellation pattern (any nontrivial
# realization of the constellation can be used) and the
# other a block pattern (the only important thing is the
# relative block sizes). For instance to compute the
# coefficient for the three term AP for the currently
# best known block configuration we use the following.
#
# Note that it stores the number as the variable CCoeff
#
blockConstellation([1,2,3],[28,6,28,37,59,116,116,59,37,28,6,28]);#
# we can also use the function to take an approximately
# optimal block structure and perturb it to find the
# best optimal block structure. For instance,
# if experimentally we found a good block structure for
# the constellation [0,1/3,1] was
# 1101,193,577,583,989,1434,1115,2833,3680,3681,2830,1113,1434,988,582,575,194,1098
# then we run the same function but first set the
# variable perturb to be true.
#
# There are several variables stored, namely:
# preCCoeff = coefficient before perturbing
# M = matrix set up to do the perturbation
# CCoeff = coefficient after perturbing
# blocks = block structure after perturbing
# Spatnew = the perturbed new subdivision pattern
#
perturb:=true:
blockConstellation([0,1/3,1],[1101,193,577,583,989,1434,1115,2833,3680,3681,2830,1113,1434,988,582,575,194,1098]);