print(` Version of April 18, 1997`):
print(`This is GJSAW, A Maple package`):
print(`that implements the Goulden-Jackson Cluster method,`):
print(`as extended and modified by John Noonan and Doron Zeilberger`):
print(`to deal with the enumeration of Self Avoiding Walks.`):
lprint(``):
print(`Written by John Noonan and Doron Zeilberger.`):
lprint(``):
print(`It is one of the packages accompanying Noonan and Zeilberger's `):
print(`article: "The Goulden-Jackson Cluster Method: Extensions,`):
print(`Applications, and Implementations", where the method`):
print(`is explained in detail, and then extended and generalized in`):
print(`various directions. `):
print(`See also John Noonan's paper: New upper bounds for connective`):
print(`constants for self-avoiding walks , preprint, available from`):
print( ` http://www.math.temple.edu/~noonan `):
lprint(``):
print(`The paper, and the packages, including the present one`):
print(`is available from http://www.math.temple.edu/~zeilberg/gj.html`):
lprint(``):
print(`For general help, and a list of the available functions,`):
print(` type "ezra();". For specific help type "ezra(procedure_name)" `):
lprint(``):
lprint(`Warning: t,s,C,x are global variables, do not use them!`):
 print(`Please report bugs to: zeilberg@math.temple.edu `):
 
 
ezra:=proc()
if args=NULL then
 
print(`This is GJSAW`):
print(`A Maple package that applies the Noonan-Zeilberger `):
print(`implementation of the Goulden-Jackson Cluster method`):
print(` to find series expansions for finite-memory SAWs`):
print(`For help, and list of procedures type "ezra();" `):
lprint(``):
print(`Written by John Noonan and Doron Zeilberger`):
lprint(``):
print(` Contains Procedures `):
print(` mist , Mist , Series , SeriesRect `):
print(`The main procedure is Series`):
print(` For specific help type "ezra(procedure_name)" `):
fi:
 
 
 
if nops([args])=1 and op(1,[args])=`SeriesRect` then
print(`SeriesRect(L): Given L, computes the sequence a(n):=`):
print(`the number of n-step walks in the plane that never traverse a`):
print(` rectangle up to L terms`):
fi:
 
 
if nops([args])=1 and op(1,[args])=`Series` then
print(`Series(MEMO,DIM,NUTERMS): Given the memory, MEMO (a positive even `):
print(`integer), the dimension, DIM, and the number of terms NUTERMS,`):
print(`finds the list whose (n+1)^(th) member is the number of memory-MEMO`):
print(`(pseudo-)SAWs in the DIM-dimensional hypercubic lattice .`):
print(`NUTERMS is the desired number of terms in the list`):
fi:
 
 
if nops([args])=1 and op(1,[args])=`GJ1` then
print(`GJ1(DIM,MISTAKES,L) finds the first L terms `):
print(`of the sequence: number of words of length 2*DIM,  in the alphabet`):
print(`\pmDIM, \pm(DIM-1), \pm1 , that avoid`):
print(`as factors,`):
print(`( i.e. subsequences of consecutive letters) the members of the `):
print(`set of lists MISTAKES and its images.`):
print(` For example to find the sequence a_n (0<=n<=15)`):
print(`where a_n:= wt. numbers of words in {-1,-2,1,2,3} without factors in `):
print(` MISTAKES={[1,1],[2,2],[-1,-1],[-2,-2]}. type GJ1(2,{[1,1]},15)`):
fi:
 
 
if nops([args])=1 and op(1,[args])=`Lclust1` then
 print(`Lclust1(MISTAKES,NUTERMS,DIM): computes the first NUTERMS terms of`):
 print(`the L-cluster expansion of the g.f. if the set of mistakes is`):
 print(` MISTAKES where MISTAKES is a set of canonical representatives`):
 print(`and the actual set of mistakes is the set of all images under`):
 print(` all signed permutations. For example Lclust1({[1,2]},10,2) `):
fi:
 
 
if nops([args])=1 and op(1,[args])=`bdok` then
 print(`bdok(MISTAKES) checks that the set of Mistakes is minimal with `):
 print(`respect to containment`):
fi:
 
 
if nops([args])=1 and op(1,[args])=`mist` then
 
print(`mist(memo,dim) : the set of canonical representatives of the set of`):
print(`mistakes with memory memo and dimension dim`):
 
fi:
 
if nops([args])=1 and op(1,[args])=`Mist` then
 
print(`Mist(memo,dim) : the set of all `):
print(`mistakes with memory memo and dimension dim`):
 
fi:
 
 
end:
 
 
 
 
#mist(memo,dim) : the set of canonical representatives of the set of
#mistakes with memory memo and dimension dim
 
squarelist:=proc(N,x) local i,b,c,L,L1:
L:=[]:
 
if N=0 then RETURN([]):
fi:
L1:=squarelist(N-1,x):
L:=[]:
for b from 1 to N do
 for c from 0 to N-1 do
L:=[op(L),[seq(x[1],i=1..N-c),seq(x[2],i=1..b),seq(-x[1],i=1..N),seq(-x[2],i=1..
b),seq(x[1],i=1..c)]]:
 od:
 for c from 0 to b-1 do
L:=[op(L),[seq(x[1],i=1..b-c),seq(x[2],i=1..N),seq(-x[1],i=1..b),seq(-x[2],i=1..
N),seq(x[1],i=1..c)]]:
 od:
od:
 
 
[op(L1),op(L)]:
end:
 
 
 
#heres the main program.mist generates all canonical mistakes up to 
#memory=mem in dimension=dim
 
mist:=proc(memo,dim) local mem,i,L:
 
if not type(memo/2,integer) or memo=0 then
 ERROR(`The first argument, memo, must be a positive even integer`):
fi:
mem:=memo/2:
L:=badlist1(mem,x,dim):
for i from 1 to dim do
  L:=subs(x[i]=i,L):
od:
L:
end:
 
 
badlist1:=proc(N,x,d) 
local ww,i,j,k,l,n,SUB,L,L1,L2,L3,LL,w,m,chk,dummy,chk2,slow:
# we make a list of all sins (up to symmetry) for a self avoiding walk with
# memory=2N. x is a variable you choose so that x[1] is the first 
# direction, 
# x[2] is the second direction etc. Output is a list of sins. the sin 
# [x[1],-x[1]] corresponds to taking one step in any direction and then 
# taking one step back.
if N=1 then RETURN([[x[1],-x[1]]]):
fi:
n:=2*N:
SUB:={x[1]=1}:
for i from 2 to N do:
SUB:={op(SUB),x[i]=1}:
od:
L:=badlist1(N-1,x,d):
L1:=[[[x[1]],2]]:
LL:=[]:
for i from 3 to n do
 L2:=[]:
 for j from 1 to nops(L1) do
  w:=op(j,L1):
  L3:=grow1(w,N,x,d):
  m:=nops(L3):
  for k from 1 to m do
   chk:=op(1,op(k,L3)):
   slow:=0:
   chk2:=convert(chk,`+`)+ww:
   for l from 1 to nops(chk2) do:
    if op(l,chk2)<>ww then slow:=slow+abs(subs(SUB,op(l,chk2))): fi:
   od:
   if slow=n-nops(chk) then
    LL:=[op(LL),op(pad(chk))]:
   else
    dummy:=0:
    for l from 1 to nops(chk)/2 do:
     if convert([op(nops(chk)-2*l+1..nops(chk),chk)],`+`)=0 then
      dummy:=1:
      break:
     fi:
    od:
    if dummy=0 then L2:=[op(L2),op(k,L3)]:
    fi:
   fi:
  od:
 od:
L1:=L2:
od:
 
[op(L),op(LL)]:
 
end:
 
 
 
grow1:=proc(w,n,x,d) local walk,cur,pre,i,L,chunk:
walk:=op(1,w):
cur:=op(2,w):
pre:=op(nops(walk),walk):
L:=[]:
if cur<=n and cur<=d then
 chunk:=[op(walk),x[cur]]:
 if checkperm(chunk)=1 then
  L:=[[chunk,cur+1]]:
 fi:
fi:
for i from 1 to cur-1 do
 if x[i]=pre or -x[i]=pre then
  chunk:=[op(walk),pre]:
  if checkperm(chunk)=1 then
   L:=[op(L),[chunk,cur]]:
  fi:
 else
  chunk:=[op(walk),x[i]]:
  if checkperm(chunk)=1 then
   L:=[op(L),[chunk,cur]]:
  fi:
  chunk:=[op(walk),-x[i]]:
  if checkperm(chunk)=1 then
   L:=[op(L),[chunk,cur]]:
  fi:
 fi:
od:
L:
end:
 
 
checkperm:=proc(tem) local k,nn,j:
 nn:=nops(tem):
 for j from 1 to nn do
  for k from j+2 to nn do
   if convert([op(j..k,tem)],`+`)=0 and k-j<>nn-1 and k-j>0 then RETURN(0): fi:
  od:
 od:
1:
end:
 
 
with(combinat):
 
 
#canon(mila):Given a word in an alphabet of integers, in terms of a list
#finds an image under the action of the symmetric group that
#is least lexicographically
 
canon:=proc(mila)
local i,gu,mu,lu,n:
lu:=[1]:
mu:=op(1,mila):
gu[mu]:=1:
gu[-mu]:=-1:
n:=2:
for i from 2 to nops(mila) do
 mu:=op(i,mila):
 if not type(gu[mu],integer) then 
  gu[mu]:=n: gu[-mu]:=-n: n:=n+1:
 fi:
 lu:=[op(lu),gu[mu]]:
od:
lu:
end:
 
 
pad:=proc(walk) local i,s,S,T,T2,check,j,dummy,l,SUB, m, n, test, zz:
SUB:={}:
m:=nops(walk):
for i from 1 to m do
 if nops(op(i,walk))=1 then 
  SUB:={op(SUB),op(i,walk)=1}:
 fi:
od:
s:=convert(zz+convert(walk,`+`),set) minus {zz};
T:=[]:
for i from 1 to nops(s) do
 check:=op(i,s):
 if nops(check)=1 then T:=[op(T),-check]:
 else
   if op(1,check)<0 then 
    for j from 1 to -op(1,check) do
     T:=[op(T),op(2,check)]:
    od:
   else
    for j from 1 to op(1,check) do
     T:=[op(T),-op(2,check)]:
    od:
   fi:
 fi:
od:
T2:=permute(T):
#print(T2,T):
n:=m+nops(T):
S:=[]:
for i from 1 to nops(T2) do
 test:=[op(walk),op(op(i,T2))]:
 dummy:=0:
 for j from 0 to nops(T)-1 do
    for l from 1 to (n-j)/2 do:
#print([op(n-j-2*l+1..n-j,test)]);
     if j<>0 or l<>n/2 then
     if convert([op(n-j-2*l+1..n-j,test)],`+`)=0 then
      dummy:=1:
      break:
     fi:
     fi:
    od:
 od:
 if dummy=0 then S:=[op(S),test]:
 fi:
od:
S:
end:
 
 
 
 
sqmist:=proc(n) local L:
L:=squarelist(n,x):
L:=subs({x[1]=1,x[2]=2},L):
end:
 
 
 
grow1:=proc(w,n,x,d) local walk,cur,pre,i,L,chunk:
walk:=op(1,w):
cur:=op(2,w):
pre:=op(nops(walk),walk):
L:=[]:
if cur<=n and cur<=d then
 chunk:=[op(walk),x[cur]]:
 if checkperm(chunk)=1 then
  L:=[[chunk,cur+1]]:
 fi:
fi:
for i from 1 to cur-1 do
 if x[i]=pre or -x[i]=pre then
  chunk:=[op(walk),pre]:
  if checkperm(chunk)=1 then
   L:=[op(L),[chunk,cur]]:
  fi:
 else
  chunk:=[op(walk),x[i]]:
  if checkperm(chunk)=1 then
   L:=[op(L),[chunk,cur]]:
  fi:
  chunk:=[op(walk),-x[i]]:
  if checkperm(chunk)=1 then
   L:=[op(L),[chunk,cur]]:
  fi:
 fi:
od:
L:
end:
 
 
 
numsuff:=proc(n) local i,L,M,nn,try1,j,lis:
 
M:=mist(n):
nn:=nops(M):
 
L:={}:
for i from 1 to nn do
 try1:=op(i,M):
 lis:=[]:
 for j from 1 to nops(try1)-1 do
  lis:=[op(lis),[op(1..j,try1)]]:
 od:
 L:=L union convert(lis,set):
od:
nops(L):
end:
 
 
#qdir(word) computes the number of dimensions that word takes on.
#thus qdir([1,2,-1,-2]) will return 2 while qdir{[1,1,1,1,2,-3]) will
# return 3.
qdir:=proc(word) local i,q,sum,n:
n:=nops(word):
sum:=1:
for i from 1 to n do
 sum:=sum+q^abs(op(i,word)):
od:
nops(sum)-1:
end:
 
 
sap:=proc(n,d) local i,M,nn,try1,ntry,L:
M:=mist(n/2,d):
nn:=nops(M):
L:={}:
for i from 1 to nn do:
 try1:=op(i,M):
 ntry:=nops(try1):
 if ntry=n then 
#  dir:=qdir(op(i,M)):
  L:={op(L),op(khaverim(try1,d))}:
#  for j from 2 to n/2 do
#   L:={op(L),canon([op(n-j+1..n,try1),op(1..n-j,try1)])}:
#  od:
#  tot:=tot+2*dir*binomial(d,dir)/n:
 fi:
od:
#print(L);
nops(L)/n/2:
end:
 
 
 
#khaverim(mila,n):Given a word mila finds the set of all its 
#images under the action of the group of signed permutations
 
khaverim:=proc(mila1,n)
local i,gu,r,lu,mila,mila2,gu1,j:
 
mila:=canon(mila1):
 
mila2:=[]:
 
for i from 1 to nops(mila) do
 mila2:=[op(mila2),abs(op(i,mila))]:
od:
 
r:=nops(convert(mila2,set)):
if n<r then
 ERROR(`Not enough letters`):
fi:
 
 
gu:=permute(n,r):
gu1:=simanim(r):
 
lu:={}:
 
for i from 1 to nops(gu) do
for j from 1 to nops(gu1) do
lu:=lu union {peula(mila,op(i,gu),op(j,gu1))}:
od:
od:
 
lu:
 
end:
 
 
#peula(mila,perm): Given a word in the alphabet {1,2,..,r} applies
#to it the (partial) permutation perm i.e. i is replaced by perm[i]
 
peula:=proc(mila,perm,siman)
local gu,i:
gu:=[]:
 
if nops(perm)<>nops(siman) then
ERROR(`Too bad`):
fi:
 
for i from 1 to nops(mila) do
gu:=[op(gu), sign(op(i,mila))*
op(abs(op(i,mila)),siman)*op(abs(op(i,mila)),perm)]:
od:
 
gu:
end:
 
 
 
simanim:=proc(r)
local i,mu,gu,evar:
 
if r=0 then
 RETURN([[]]):
fi:
 
gu:=[]:
 
mu:=simanim(r-1):
 
for i from 1 to nops(mu) do
 evar:=op(i,mu):
gu:=[op(gu),[op(evar),1],[op(evar),-1]]:
od:
gu:
end:
 
 
Mist:=proc(memo,dim) 
local gu,gu1,i:
 
gu:=mist(memo,dim) :
gu1:={}:
 
for i from 1 to nops(gu) do
gu1:=gu1 union khaverim(op(i,gu),dim):
od:
 
gu1:
 
end:
 
 
 
 
 
#bdok(MISTAKES) checks that the set of Mistakes is minimal with 
#respect to containment
 
bdok:=proc(MISTAKES)
local mila,i,j,k,mila1:
 
for i from 1 to nops(MISTAKES) do
mila:=op(i,MISTAKES):
 
for j from 1 to nops(mila) do
 for k from j to nops(mila) do
  mila1:=[op(j..k,mila)]:
  if k-j<nops(mila)-1 and member(mila1,MISTAKES) then
   RETURN(0):
  fi:
 od:
od:
od:
1:
end:
 
 
 
 
hafukh:=proc(mila)
local i,mila1:
mila1:=[]:
 
for i from   nops(mila) by -1  to 1  do
mila1:=[op(mila1),op(i,mila)]:
od:
mila1:
end:
 
 
 
canon1:=proc(mila)
hafukh(canon(hafukh(mila))): 
end:
 
 
 
Lclust1:=proc(MISTAKES1,NUTERMS,DIM)
local LEN,SUFFI,i,j,x,Lt,LIS,mila,I1,r,pref,tot,gu,resh,i1,suffi,
MISTAKES,kapi1,kapi2,mila1:
 
MISTAKES:={}:
for i from 1 to nops(MISTAKES1) do
 mila:=canon1(op(i,MISTAKES1)):
kapi1[op(mila)]:=kama(mila,DIM):
MISTAKES:=MISTAKES union {mila}:
od:
 
 
LIS:=[]:
 
SUFFI:={[]}:
 
for i from 1 to nops(MISTAKES) do
 mila:=op(i,MISTAKES):
 
 for j from 2 to nops(mila) do
 mila1:=[op(j..nops(mila),mila)]:
kapi2[op(mila1)]:=kama(mila1,DIM):
  SUFFI:=SUFFI union {mila1}:
 od:
od:
 
LEN:=0:
 
for i from 1 to nops(MISTAKES) do
  mila:=op(i,MISTAKES):
  if nops(mila)>LEN then
    LEN:=nops(mila):
  fi:
od:
 
 
for i from 1 to nops(SUFFI) do
mila:=op(i,SUFFI):
Lt[op(mila)]:=[seq(0,j=1..LEN-nops(mila))]:
od:
 
 
for I1 from 0 to NUTERMS do
 
tot:=0:
 
for i from 1 to nops(SUFFI) do
mila:=op(i,SUFFI):
x[op(mila)]:=0:
od:
 
for i from 1 to nops(MISTAKES) do
 
  mila:=op(i,MISTAKES):
 
  gu:=0:
 
  if nops(mila)=I1 then
   gu:=-1:
  fi:
 
  for r from 1 to nops(mila)-1 do
    pref:=canon1([op(1..nops(mila)-r,mila)]):
     
     if member(pref,SUFFI) then
      gu:=gu-op(nops(Lt[op(pref)])-r+1,Lt[op(pref)]):
     fi:
 
  od:
 
 
  for j from 2 to nops(mila) do
   suffi:=[op(j..nops(mila),mila)]:
   x[op(suffi)]:=x[op(suffi)]+kapi1[op(mila)]
/kapi2[op(suffi)]*gu:
  od:
 
 
tot:=tot+kapi1[op(mila)]*gu:
 
od:
 
 
 
LIS:=[op(LIS),tot]:
 
for i1 from 1 to nops(SUFFI) do
mila:=op(i1,SUFFI):
resh:=Lt[op(mila)]:
Lt[op(mila)]:=[op(2..nops(resh),resh),x[op(mila)]]:
od:
 
od:
 
LIS:
 
end:
 
 
GJ1:=proc(DIM,MISTAKES,L)
local NUMLETTERS,lu,resh,i,gu,t,resh1:
 
NUMLETTERS:=2*DIM:
 
resh:=Lclust1(MISTAKES,L+2,DIM):
 
lu:=0:
 
for i from 0 to L+1 do
lu:=lu+op(i+1,resh)*t^i:
od:
 
 
gu:=1/(1-NUMLETTERS*t-lu):
 
gu:=taylor(gu,t=0,L+1):
 
resh1:=[]:
 
for i from 0 to L do
resh1:=[op(resh1),coeff(gu,t,i)]:
od:
 
resh1:
end:
 
 
 
#Series(MEMO,DIM,NUTERMS): Given the memory, MEMO (a positive even 
#integer), the dimension, DIM, and the number of terms NUTERMS,
#finds the list whose (n+1)^(th) member is the number of memory-MEMO
#(pseudo-)SAWs in the DIM-dimensional hypercubic lattice .
#NUTERMS is the desired number of terms in the list
 
Series:=proc(MEMO,DIM,NUTERMS)
local MISTAKES:
MISTAKES:=mist(MEMO,DIM) :
GJ1(DIM,MISTAKES,NUTERMS):
end:
 
 
kama:=proc(mila,n)
local i,r,mila2:
 
 
mila2:=[]:
 
for i from 1 to nops(mila) do
 mila2:=[op(mila2),abs(op(i,mila))]:
od:
 
r:=nops(convert(mila2,set)):
if n<r then
 ERROR(`Not enough letters`):
fi:
 
2^r*n!/(n-r)!:
 
end:
 
 
 
 
#rectmis(L) all represenatives of rectangular mistakes
 
rectmis:=proc(L)
local gu,a,b,i,mila,i1:
gu:={[1,-1]}:
 
for a from 1 to trunc(L/2)  do
for b from 1 to trunc((L-2*a)/2) do
for i from 1 to a do
 mila:=[seq(1,i1=1..i),seq(2,i1=1..b),seq(-1,i1=1..a),seq(-2,i1=1..b),
seq(1,i1=1..a-i)]:
gu:=gu union {mila}:
od:
od:
od:
gu:
end:
 
 
 
#SeriesRect(L): Given L, computes the sequence a(n):=
#the number of n-step walks in the plane that never traverse a rectangle
#up to L terms
 
SeriesRect:=proc(L)
local MISTAKES:
MISTAKES:=rectmis(L+1):
GJ1(2,MISTAKES,L):
end: