Find Group (Up to size 200)# Input: TDSteps (Two Dimensional Step Set)
# Output: List (List of Group Elements)
genGP := proc(TDSteps)
local add2,Am1,A0,A1,Bm1,B0,B1,N,S,phi,psi,List,t2I,t2P,newI,newP,i,k,tst;
# Various things needed later
add2 := set -> add(X^k[1]*Y^k[2],k=set):
S := add2(TDSteps):
Am1 := coeff(S,Y,-1): A0 := coeff(S,Y,0): A1 := coeff(S,Y,1):
Bm1 := coeff(S,X,-1): B0 := coeff(S,X,0): B1 := coeff(S,X,1):
phi := S -> simplify(subs(X=S[1],Y=S[2],[(1/X)*Bm1/B1,Y])):
psi := S -> simplify(subs(X=S[1],Y=S[2], [X,(1/Y)*Am1/A1])):
# Assume we know the group is finite, then eventually this iteration will stabalize
List:= [ [[],[x,y]] ]:
N := -1:
while (N <> nops(List)) and (nops(List)<200) do
N := nops(List):
for i in List do
newI := true:
newP := true:
t2I := phi(i[2]):
t2P := psi(i[2]):
t2I := map(x->normal(expand(numer(x))/expand(denom(x))),t2I):
t2P := map(x->normal(expand(numer(x))/expand(denom(x))),t2P):
for k in List do
if k[2]=t2I then newI := false: fi:
if k[2]=t2P then newP := false: fi:
if not (newI or newP) then break: fi:
od:
if newI then List := [ [[phi,op(i[1])],t2I], op(List)]: fi;
if newP then List := [ [[psi,op(i[1])],t2P], op(List)]: fi;
if nops(List)>200 then break: fi:
od:
od:
if
nops(List) < 200 then return List;
else
return FAIL;
fi;
end:Find Orbit Sum# Input: TDSteps (Two Dimensional Step Set)
# Output: Ex (Expression resulting from orbit sum)
# Requires GenGP2D function defined above
OrbitSum2D := proc(TDSteps) local List,EQ,Ex,T,K,LST,i,k;
K := 1-simplify(t*add(x^k[1]*y^k[2],k=TDSteps)):
List := genGP(TDSteps):
# Sum the Kernel Equation after it's acted upon by the group elements
EQ := Kk*X*Y*Q(X,Y) = X*Y + F(Y) + G(X):
Ex := 0:
for i from 1 to nops(List) do
T := subs(X=List[i][2][1],Y=List[i][2][2],EQ):
Ex := Ex+(-1)^nops(List[i][1])*T:
od:
#The result:
Ex := expand(lhs(Ex)/(Kk*x*y)) = simplify(rhs(%)/(K*x*y)):
return Ex;
end:Code to Generate 2D Counting Sequences# Updates walk locations after one iteration
# Input: Current Locs, Place to Store New Locs, Step Set, and Bounds
Update2D := proc(Locations,NewLocations,Steps,max) local i,j,s;
for i from 0 to max do
for j from 0 to max do
if type(Locations[i,j],integer) then
for s in Steps do
if (i+s[1])<0 or (j+s[2])<0 then ;
elif type(NewLocations[i+s[1],j+s[2]],integer) then
NewLocations[i+s[1],j+s[2]] := NewLocations[i+s[1],j+s[2]]+Locations[i,j];
else
NewLocations[i+s[1],j+s[2]] := Locations[i,j] #Was previously zero
end if;
end do;
end if;
end do;
end do;
end:
# Sums terms in a table
# Input: Table and Bounds
SumTab2D := proc(Table,max) local i,j,out;
out := 0;
for i from 0 to max do
for j from 0 to max do
if type(Table[i,j],integer) then out := out + Table[i,j];
fi;od;od;
return out;
end:
# Returns list with the number of walks up to length k
# Input: Step set (as list of lists) and number of steps
Count2D := proc(StepsT,k) local Steps,tt,Locations,NewLocations,F,Points,temp;
Steps := eval(StepsT):
Locations := 'Locations':
Locations[0,0] := 1:
NewLocations := 'NewLocations':
F := [1]:
for temp from 1 to k do
Update2D(Locations,NewLocations,Steps,temp-1);
tt := SumTab2D(NewLocations,temp):
F := [op(F),tt];
Locations := copy(NewLocations):
NewLocations := 'NewLocations':
end do:
return F:
end:Code to Generate 3D Counting Sequences# Updates walk locations after one iteration
# Input: Current Locs, Place to Store New Locs, Step Set, and Bounds
Update3D := proc(Locations,NewLocations,Steps,max) local i,j,k,s;
for i from 0 to max do
for j from 0 to max do
for k from 0 to max do
if type(Locations[i,j,k],integer) then
for s in Steps do
if (i+s[1])<0 or (j+s[2])<0 or (k+s[3])<0 then ;
elif type(NewLocations[i+s[1],j+s[2],k+s[3]],integer) then
NewLocations[i+s[1],j+s[2],k+s[3]] := NewLocations[i+s[1],j+s[2],k+s[3]]+Locations[i,j,k];
else
NewLocations[i+s[1],j+s[2],k+s[3]] := Locations[i,j,k] #Was previously zero
end if;
end do;
end if;
end do;
end do;
end do;
end:
# Sums terms in a table
# Input: Table and Bounds
SumTab3D := proc(Table,max) local i,j,k,out;
out := 0;
for i from 0 to max do
for j from 0 to max do
for k from 0 to max do
if type(Table[i,j,k],integer) then out := out + Table[i,j,k];
fi;od;od;od;
return out;
end:
# Returns list with the number of walks up to length k
# Input: Step set (as list of lists) and number of steps
Count3D := proc(StepsT,k) local Steps,tt,Locations,NewLocations,F,Points,temp;
Steps := eval(StepsT):
Locations := 'Locations':
Locations[0,0,0] := 1:
NewLocations := 'NewLocations':
F := [1]:
for temp from 1 to k do
Update3D(Locations,NewLocations,Steps,temp-1);
tt := SumTab3D(NewLocations,temp):
F := [op(F),tt];
Locations := copy(NewLocations):
NewLocations := 'NewLocations':
end do:
return F:
end:read "TwoDSteps.txt":# The list Reducible holds all models which are two dimensional
nops(Reducible);# The step sets are arranged so that the z coordinate is redundant:
Reducible[1234];
Count3D(%,15);
map(x->x[1..2],Reducible[1234]);
Count2D(%,15);# The list TwoDSteps holds all the unique 2D step sets which arise
nops(TwoDSteps);# The list ProvenOrbit2D holds the 3D models with non-zero orbit sum which are
# proven to admit D-finite generating functions by the kernel method
nops(ProvenOrbit2D);# ModelS1 is the step set with non-zero orbit sum which has an error in coefficient extraction.
# It is proven to have a D-finite generating function through a Computer Algebra approach
ModelS1;
rhs(OrbitSum2D(ModelS1));# ModelS1bar is the step set with zero orbit sum proven
# algebraic by a Computer Algebra approach
ModelS1bar;
rhs(OrbitSum2D(ModelS1bar));# The list ZeroOrbit holds the remaining models with zero orbit sum
nops(ZeroOrbit);
OrbitSum2D(ZeroOrbit[1]);# The following code can calculate the positive part of a rational function
# This is very temperamental in Maple.
# It is suggested to use increase the accuracy N if undesired behavior occurs
N := 10:
EX := rhs(OrbitSum2D(ProvenOrbit2D[12]));
temp := series(EX,t,2*N+1):
temp := map(a -> series(a,x,2*N+1), convert(temp,polynom)):
temp := map(a -> series(a,y,2*N+1), convert(temp,polynom)):
temp := convert(temp,polynom):
temp := add(add(add(x^j*y^k*t^i*coeff(coeff(coeff(temp,t,i),x,j),y,k),i=0..N+1),j=0..N+1),k=0..N+1):
temp := subs(x=1,y=1,temp):
[seq(coeff(temp,t,k),k=0..N)];# This matches the counting sequence of the model, as desired
Count2D(ProvenOrbit2D[12],10);LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYjLUkjbWlHRiQ2I1EhRic=