print(` Version of OCT 6 2003`):
print(`This is SymbolicGJ, A Maple package`):
print(`that symbolicly implements the Goulden-Jackson Cluster method .`):
print(`Written by Xiangdong Wen and Doron Zeilberger.`):
print(`part of the functions are copied from:`): 
print(`http://www.math.rutgers.edu/~zeilberg/gj.html`):
lprint(``):
print(`Ian Goulden and David Jackson's original treatment can be found `):
print(`in their paper that appeared in J. London Math. Soc.(2)20(1979),`):
print(` 567-576. It is also described in their book`):
print(`Combinatorial Enumeration, John  Wiley, 1983, New York`):
print(`See also R. Stanley's Enumerative Combinatorics, v.1, 266, 280.`):
print(`and Guibas, Odlyzko,JCT(A)30(1981),183-208.`):
lprint(``):
print(`For general help, and a list of the available functions,`):
print(` type ezra();. For specific help type ezra(procedure_name) `):
lprint(``):

ezra:=proc()
if args=NULL then
print(`SYMBOLIC_GJ`):
print(`A Maple package that Symbolicly implements`):
print(` the Goulden-Jacskon Cluster method`):
print(` for finding generating functions and series expansions`):
print(`for the number of words avoiding a prescribed, finite, set of`):
print(`mistakes. Contains Procedures: `):
print(` SGJ,SSGJ,SAW_SSGJ`):
print(` For specific help type ezra(procedure_name) `):
fi:
   
if nops([args])=1 and op(1,[args])=`SGJ` then
print(`SGJ(n,MISTAKES,s) finds the generating function`):
print(`g(n,s), is the g.f. sum(a[i]*s^i,i=0..infty)  `):
print(`where a[i]:=number of words of length i without any mistakes`):
print(`from MISTAKES, and n is a symbol indicate the total number of letters.`):
print(`For example to get the g.f. for a_m:=number of words of length m`):
print(`in the alphabet {1,2,3,...,n}, avoiding 123,231,312, type`):
print(` SGJ(n,{[1,2,3],[2,3,1],[3,1,2]},s)`):
fi:

if nops([args])=1 and op(1,[args])=`SSGJ` then
print(`SSGJ(n,MISTAKES,s) finds the generating function`):
print(`g(s), is the g.f. sum(a[i]*s^i,i=0..infty)  `):
print(`where a[i]:=number of words of length i without any mistakes`):
print(`from MISTAKES and its images.`):
print(`where n is a symbol indicates the total number of letters`):
print(`For example to get the g.f. for a_m:=number of words of length m`):
print(`in the alphabet {1,2,3,...,n}, avoiding 123,132,231,213,312, and 321 type`):
print(` SSGJ(n,[[1,2,3]],s);`):
fi:

if nops([args])=1 and op(1,[args])=`SAW_SSGJ` then
print(`SAW_SSGJ(n,MISTAKES,s)`):
print(`Finding generating functions `):
print(`for the number of words, in the alphabet -1,1,-2,2,..,-d,d`):
print(` avoiding a prescribed, finite, set of "MISTAKES"`):
print(`the set of mistakes is closed under the action of the group of`):
print(`signed permutations`):
print(`For example to get the g.f. for a_m:=number of words of length m`):
print(`in the alphabet {-1,1,-2,2,-3,3,...,-n,n}, 2 memory self avoiding walk,type `):
print(` SAW_SSGJ(n,[[1,-1]],s);`):

fi:
end:

########################################################################################
########################################################################################
##
##This part is for SSGJ
##
########################################################################################
########################################################################################

######################################################################################## 
###isostring(u,v),judges whether u and v are isostrings,
###if u and v are iso return true, otherwise false.
###usage: isostring([1,2,3,4,5],[5,3,4,2,1]);

isostring:=proc(u,v)
local i,j,t:
if(nops(u)<>nops(v))then
return false;
fi:

for i from 1 to nops(u) do
t:={op(i,v)};
for j from i+1 to nops(u)do
if op(i,u)=op(j,u) then
t:=t union {op(j,v)};
fi:
od:
if nops(t)>=2 then return false;
fi:
od:

for i from 1 to nops(v) do
t:={op(i,u)};
for j from i+1 to nops(v)do
if op(i,v)=op(j,v) then
t:=t union {op(j,u)};
fi:
od:
if nops(t)>=2 then return false;
fi:
od:
return true:
end:

##########################################################################################
###findequ(n,u,seqv,s) get the equations for string u with strings seqv.
###usage: findequ(n,[1,2,3],[[1,2,3],[1,1]],s)

findequ:=proc(n,u,seqv,s)
local i,j,h,t,x,y,m,term,equ,k,v,vset,yset:
equ:=-s^(nops(u)):

for h from 1 to nops(seqv) do
v:=op(h,seqv):

for i from 1 to nops(u) do
x:=[op(i..nops(u),u)]:
if nops(v)> nops(u)-i+1 then
y:=[op(1..nops(u)-i+1,v)]:

if isostring(x,y) then

vset:=convert(v,set);
yset:=convert(y,set); 
m:=nops(convert(v,set) minus yset);
k:=nops(yset);
term:=C(op(v)):

for j from 1 to m do
term:=term*(n-k-j+1): 
od:
equ:=equ-term*s^(nops(u)-nops(y)):
fi:
fi:
od:
od:
C(op(u))=equ:
end:

#####################################################################################
###kernalarr(Arr) clean Arr with only one copy of iso string left for each type.
###usage:kernalarr([[1,2,3],[3,2,1],[1,1]]);

kernalarr:=proc(Arr)
local i,j,ans,isin:
ans:=[];
for i from 1 to nops(Arr) do
  isin:=0:
  for j from 1 to nops(ans) do
    if isostring(op(j,ans),op(i,Arr)) then
       isin:=1;
    fi:
  od:
  if isin=0 then
  ans:=[op(ans),op(i,Arr)]:
  fi:
od:
ans:
end:

#############################################################################################
###IsIsoSubString(arr1,arr2) test whether arr2 is a sub string of arr1.
###usage:IsIsoSubString([1,2,3,4],[4,3,2]);

IsIsoSubString:=proc(arr1,arr2)
local i,j;
if nops(arr2)>nops(arr1) then
return false:
fi:

for i from 1 to nops(arr1)-nops(arr2)+1 do
   if isostring(arr2,[op(i..i+nops(arr2)-1,arr1)]) then
   return true:
   fi:
od:
return false:
end:
##############################################################################################
###GetOffIsoSupstring(arr) clean the string arr.
###usage:GetOffIsoSupstring([[1,2,2,3],[2,2]]);

GetOffIsoSupstring:=proc(arr)
local i,j,k,isbigstring,ans;
ans:=[]:
for i from 1 to nops(arr) do
   isbigstring:=0:
   for j from 1 to nops(arr) do

   if i<>j then   
   if IsIsoSubString(op(i,arr),op(j,arr)) then
   isbigstring:=1:
   fi:
   fi:
   od:
   if isbigstring=0 then
   ans:=[op(ans),op(i,arr)]:
   fi:
od:
ans:
end:

##########################################################################################
####SSGJ(n,mistake,s), symbolic symetric goulden Jackson, where n and s are symboles
####mistake are arrays of mistakes,
####usage:SSGJ(n,[[1,2,3],[1,1]],s);

SSGJ:=proc(n,mistake,s)
local i,j,k,eq,var,res,lu,v,zm,coef,ans,mistakes:

mistakes:=kernalarr(mistake);
mistakes:=GetOffIsoSupstring(mistakes):
eq:={}:
var:={}:
for i from 1 to nops(mistakes) do
eq:= eq union {findequ(n,op(i,mistakes),mistakes,s)}:
var:= var union {C(op(op(i,mistakes)))}:
od:

var:=solve(eq,var):

lu:=1-s*n:

for i from 1 to nops(mistakes) do
 v:=op(i,mistakes):
  zm:=nops(convert(v,set)):

  coef:=1:
  for j from 1 to zm do
  coef:=coef*(n-j+1):
  od:

 lu:=lu-coef*subs(var,C(op(v))):
od:
 
ans:=normal(1/lu):

collect(numer(ans),s)/collect(denom(ans),s) :
end:


########################################################################################
########################################################################################
##
##This part is for SGJ
##
########################################################################################
########################################################################################

########################################################################################
###SGJ, symbolic goulden-jackson cluster method.
###usage:SGJ(n,{[1,2,3],[1,1]},s);

SGJ:=proc(n,MISTAKES1,s)
local v,eq, var,i,lu,C,MISTAKES,ans:
 
MISTAKES:=Hakten(MISTAKES1):
eq:={}:
var:={}:
 
for i from 1 to nops(MISTAKES) do
 v:=op(i,MISTAKES):
 eq:= eq union {findeqz(v,MISTAKES,C,s)}:
 var:=var union {C[op(v)]}:
od:
 
var:=solve(eq,var):
 
 lu:=1-n*s:
 
for i from 1 to nops(MISTAKES) do
 v:=op(i,MISTAKES):
 lu:=lu-subs(var,C[op(v)]):
od:
 
ans:=normal(1/lu):

collect(numer(ans),s)/collect(denom(ans),s) :

end:

################################################################################
###findeqz sets up the equ C[v]= s+t*Sum_u overlap(u,v,x) *C[u]
###usage: findeqz([1,2,3],{[1,2,3],[1,1,1]},C,s);

findeqz:=proc(v,MISTAKES,C,s)
local eq,i,u:
eq:=-1:
 
for i from 1 to nops(v) do
 eq:=eq*s:
od:
 
for i from 1 to nops(MISTAKES) do
 u:=op(i,MISTAKES):
 eq:=eq-overlapz(u,v,s)*C[op(u)]:
od:
 
C[op(v)]-eq=0:
 
end:
###############################################################################
###overlapz is a procedure that given two words u and v, and a variable s
###computes the weight-enumerator of all v\suffix(u),
###for all suffixes of u that are prefixes of v, but with uniform weight s
###usage: overlapz([1,2,3,4],[3,4,5,6],s);

overlapz:=proc(u,v,s)
local i,j,lu,gug:
 
lu:=0:
 
for i from 2 to nops(u) do
  for j from i to nops(u) while (j-i+1<=nops(v) and op(j,u)=op(j-i+1,v))
         do
   :
  od:   
 
if j-i=nops(v) and u<>v then
 ERROR(v,`is a subword of`,u,`illegal input`):
 fi:
 
 
 
 if j=nops(u)+1 and (i>1 or j>2) then
  gug:=1:
 
   gug:=gug*s^(nops(v)-(j-i)):
 
  lu:=lu+gug:
 
fi:
 
od:
 
lu:
 
end:


########################################################################################## 
###Hakten1(B) removes all superflous words
###usage: hakten1({[1,2,3],[0,1,2,3,4]});

hakten1:=proc(B)
local w,i:
 
for i from 1 to nops(B) do
w:=op(i,B):
 
if superflous(B,w)=1 then
 RETURN(B minus {w}):
fi:
od:
B:
end:
########################################################################################## 
###Hakten(B) removes all superflous words
###usage: Hakten({[1,2,3],[0,1,2,3,4]});

Hakten:=proc(B)
local B1,B2:
 
B1:=B:
 
B2:=hakten1(B1):
 
while B1<>B2 do
B1:=B2:
B2:=hakten1(B2):
od:
B2:
end:


###########################################################################################
###superflous(B,w) given a set of bad words, and one of its
###members, finds whether it is superflous,
###that is, whether it is a factor of another one, and deletes the larger one
###usage:superflous({[1,2,3],[1,2,3,4]},[1,2,3]);

superflous:=proc(B,w)
local i,v:
 
if not type(B,set) or not type(w,list) then
 ERROR(`Bad input`):
fi:
 
if not member(w,B) then
 ERROR(`Second argument must belong to first arg.`):
fi:
 
for i from 1 to nops(B) do
v:=op(i,B):
if v<>w and issubword(w,v)=1 then
 RETURN(1):
fi:
od:
0:
end:

########################################################################################## 
##issubword(u,v) returns 1 if v is a subword of u, otherwise 0
##usage: issubword([1,2,3,4],[2,3,4]);

issubword:=proc(u,v)
local i,j:
for i from 1 to nops(u) do
  for j from i to nops(u) while (j-i+1<=nops(v) and op(j,u)=op(j-i+1,v))
         do
   :
  od:   
  if j-i=nops(v)  then
  RETURN(1):
  fi:
od:
 
0:
end:



##########################################################################################
##########################################################################################
###
###
###This part is for self avoiding walk.
###
##########################################################################################
##########################################################################################


########################################################################################### 
###saw_isostring(u,v),judges whether u and v are isostrings,
###if u and v are iso return true, otherwise false.
###usage: saw_isostring([1,2,-1,-2],[2,1,-2,-1]);

saw_isostring:=proc(u,v)
local i,j,t:
if(nops(u)<>nops(v))then
return false;
fi:

for i from 1 to nops(u) do
t:={op(i,v)};
for j from i+1 to nops(u)do
if op(i,u)=op(j,u) then
t:=t union {op(j,v)};
fi:
if op(i,u)=-op(j,u) then
t:=t union {-op(j,v)};
fi:
od:
if nops(t)>=2 then return false;
fi:
od:

for i from 1 to nops(v) do
t:={op(i,u)};
for j from i+1 to nops(v)do
if op(i,v)=op(j,v) then
t:=t union {op(j,u)};
fi:
if op(i,v)=-op(j,v) then
t:=t union {-op(j,u)};
fi:
od:
if nops(t)>=2 then return false;
fi:
od:
return true:
end:

##########################################################################################
###saw_findequ(n,u,seqv,s) get the equations for string u with strings seqv.
###usage: saw_findequ(n,[1,2,3],[[1,2,3],[1,1]],s)

saw_findequ:=proc(n,u,seqv,s)
local i,j,h,t,x,y,m,term,equ,k,v,vset,yset:
equ:=-s^(nops(u)):

for h from 1 to nops(seqv) do
v:=op(h,seqv):

for i from 1 to nops(u) do
x:=[op(i..nops(u),u)]:
if nops(v)> nops(u)-i+1 then
y:=[op(1..nops(u)-i+1,v)]:

if saw_isostring(x,y) then

#vset:=convert(v,set);
#yset:=convert(y,set); 

m:=nops(saw_kernal(v))-nops(saw_kernal(y));
k:=nops(saw_kernal(y));
term:=C(op(v)):

for j from 1 to m do
term:=term*(n-k-j+1)*2: 
od:
equ:=equ-term*s^(nops(u)-nops(y)):
fi:
fi:
od:
od:
C(op(u))=equ:
end:

#####################################################################################
###saw_kernalarr(Arr) clean Arr with only one copy of iso string left for each type.
###usage:saw_kernalarr([[1,2,3],[3,2,1],[1,1]]);

saw_kernalarr:=proc(Arr)
local i,j,ans,isin:
ans:=[];
for i from 1 to nops(Arr) do
  isin:=0:
  for j from 1 to nops(ans) do
    if saw_isostring(op(j,ans),op(i,Arr)) then
       isin:=1;
    fi:
  od:
  if isin=0 then
  ans:=[op(ans),op(i,Arr)]:
  fi:
od:
ans:
end:

#############################################################################################
###saw_kernal(Arr):

saw_kernal:=proc(Arr)
local i,j,k,ans;

ans:=[]:
for i from 1 to nops(Arr) do
   k:=0;
   for j from 1 to nops(ans) do
      if op(i,Arr)=op(j,ans) or op(i,Arr)=-op(j,ans) then
         k:=1;
      fi:
   od:
   
   if k=0 then
     ans:= [op(ans),op(i,Arr)]:
   fi:
od:
ans:
end:
#############################################################################################
###saw_IsIsoSubString(arr1,arr2) test whether arr2 is a sub string of arr1.
###usage:saw_IsIsoSubString([1,2,3,4],[4,3,2]);

saw_IsIsoSubString:=proc(arr1,arr2)
local i,j;
if nops(arr2)>nops(arr1) then
return false:
fi:

for i from 1 to nops(arr1)-nops(arr2)+1 do
   if saw_isostring(arr2,[op(i..i+nops(arr2)-1,arr1)]) then
   return true:
   fi:
od:
return false:
end:
##############################################################################################
###saw_GetOffIsoSupstring(arr) clean the string arr.
###usage:saw_GetOffIsoSupstring([[1,2,2,3],[2,2]]);

saw_GetOffIsoSupstring:=proc(arr)
local i,j,k,isbigstring,ans;
ans:=[]:
for i from 1 to nops(arr) do
   isbigstring:=0:
   for j from 1 to nops(arr) do

   if i<>j then   
   if saw_IsIsoSubString(op(i,arr),op(j,arr)) then
   isbigstring:=1:
   fi:
   fi:
   od:
   if isbigstring=0 then
   ans:=[op(ans),op(i,arr)]:
   fi:
od:
ans:
end:

##########################################################################################
####saw_SSGJ(n,mistake,s), symbolic symetric goulden Jackson, where n and s are symboles
####mistake are arrays of mistakes,
####usage:saw_SSGJ(n,[[1,-1],[1,2,-1,-2]],s);

SAW_SSGJ:=proc(n,mistake,s)
local i,j,k,eq,var,res,lu,v,zm,coef,ans,mistakes:

mistakes:=saw_kernalarr(mistake);
mistakes:=saw_GetOffIsoSupstring(mistakes):
eq:={}:
var:={}:
for i from 1 to nops(mistakes) do
eq:= eq union {saw_findequ(n,op(i,mistakes),mistakes,s)}:
var:= var union {C(op(op(i,mistakes)))}:
od:

var:=solve(eq,var):

lu:=1-2*s*n:

for i from 1 to nops(mistakes) do
 v:=op(i,mistakes):
  zm:=nops(saw_kernal(v)):

  coef:=1:
  for j from 1 to zm do
  coef:=coef*(n-j+1)*2:
  od:

 lu:=lu-coef*subs(var,C(op(v))):
od:
 
ans:=normal(1/lu):

collect(numer(ans),s)/collect(denom(ans),s) :
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: