######################################################################
##ZOO: Save this file as ZOO. To use it, stay in the                 #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read ZOO : <Enter>                                 #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Rutgers University ,                  #
#zeilberg@math.rutgers.edu.                                          # 
######################################################################
 
#Created: Nov. 2000,
#This version: Nov. 2000
#ZOO: A Maple package that counts k-board lattice animals
#using the Umbral Transfer Matrix Theorem
#It accompanies Doron Zeilberger's article:
#The Umbral Transfer Matrix Method. III. Counting Animals
#Please report bugs to zeilberg@math.rutgers.edu
 
 
print(`ZOO: A Maple package that counts k-board lattice animals`):
print(`using the Umbral Transfer Matrix Theorem`):
print(`It accompanies Doron Zeilberger's article:`):
print(`The Umbral Transfer Matrix Method. III. Counting Animals`):
lprint(``):
print(`Created: Nov. 2000`):
print(`This version: Nov. 2000`):
lprint(``):
print(`Written by Doron Zeilberger, zeilberg@math.temple.edu`):
lprint(``):
print(`Please report bugs to zeilberg@math.temple.edu`):
lprint(``):
 print(`The most current version of this  package and paper`):
 print(` are  available from`):
 print(`http://www.math.temple.edu/~zeilberg/`):
 print(`For a list of the procedures type ezra(), for help with`):
 print(`a specific procedure, type ezra(procedure_name)`):
print(`AUS1 and AUS2 are the pre-computed Umbral Schemes for`):
print(`1-board and 2-board lattice animals`):
 print(``):
 
ezra:=proc()
if args=NULL then
print(`AUS1 and AUS2 are the pre-computed Umbral Schemes for`):
print(`1-board and 2-board lattice animals`):
 print(`Contains the following procedures: Alphabet, AniUmSc, `):
  print(` ApplyUmSc, APU, `):
 print(` GapInterfaces,  Ini, Interfaces `):
 print(` IsLegalInfc , LeftLetter, LegalInterfaces, PreUmbra, `):
 print(`  PzA ,ToUmbra,`):
print(``):
 
fi:
 
if nops([args])=1 and op(1,[args])=Alphabet then
print(` Alphabet(L): The alphabet of 2D-lattice animals that have`):
print(` the property that each vertical cross-section has at most`):
print(` L boards (and hence L-1 gaps) `):
print(`For example, to get the alphabet for <=2 boards do: Alphabet(2); `)
fi:
 
if nops([args])=1 and op(1,[args])=AniUmSc then
print(`AniUmSc(m,x,y,q): An Umbral Scheme in terms of`):
print(` A 4-list whose members are `):
print(`a list of lists of Umbra, followed by the list of`):
print(`arguments of the corresponding functions`):
print(`for m-board lattice animals`):
print(`followed by the initial values`):
print(`followed by the set of terminal letters`):
fi:
 
 
if nops([args])=1 and op(1,[args])=ApplyUmSc then
print(` ApplyUmSc(UmSch,q,n,var): Given an Umbral Scheme, UmSch, finds`):
print(`the generating functions of creatures up to the wt. q^n`):
print(`followed by substition of all x and y variables to 1 and summing`):
print(`over leftmost letters`):
fi:
 
if nops([args])=1 and op(1,[args])=APU then
print(`APU(Intfc1,Gintfc1,m,n,x,y,q,a,b): Given an interface, Intfc1 for the`):
print(`m boards, and an interface Gintfc1, for the m+1 gaps,`):
print(`returns the Atomic Pre Umbra, corresponding to these`):
print(`Interfaces, i.e. the image of x[1]^a[1]*y[1]^b[1]*x[2]^a[2]...`):
print(`the right letter has n boards`):
print(` q is the variable that keeps track of the number of cells`):
print(` For example, try: `):
print(`APU([[1,2],[3,4]],[[1,1],[2,3],[4,5]],2,2,x,y,q,a,b);`):
fi:
 
if nops([args])=1 and op(1,[args])=GapInterfaces then
print(`GapInterfaces(Intfc,n): Given an interface for the boards`):
print(`in terms of a list of intervals, finds the set of`):
print(`compatible interfaces for the gaps`):
print(`Here n is the number of boards of the right letter`):
print(`For example, try GapInterfaces([[1,2],[3,4]],2);`): 
fi:
 
if nops([args])=1 and op(1,[args])=Ini then
print(`Ini(x,y,q,m): The initial weight of a leftmost letter`):
print(`with m boards according to x (for the boards), y (for the gaps)`):
print(`and q (for the length of the boards)`):
fi:
 
if nops([args])=1 and op(1,[args])=Interfaces then
print(`Interfaces(m,n): All possible interfaces between an`):
print(`n-board animal-letter from the right and an m-board`):
print(`animal letter on its left`):
print(` For example, try Interfaces(2,3); `):
fi:
 
if nops([args])=1 and op(1,[args])=IsLegalInfc then
print(`IsLegalInfc(LETTER1,interface1): given a letter LETTER1`):
print(`and an interface, interface1, that defines`):
print(`a letter to the left, returns true if the interface`):
print(`is a legal continuation`):
print(` For example, try: IsLegalInfc({{1},{2}},[[1,1],[2,3]]);`):
fi:
 
if nops([args])=1 and op(1,[args])=LeftLetter then
print(`LeftLetter(Intfc1,LETTERr): Given an interface`):
print(`Intfc1, and a letter to its right, LETTERr,`):
print(`finds the letter than results on the left`):
print(` For example, try: LeftLetter([[1,3]],{{1},{2}}); `):
fi:
 
if nops([args])=1 and op(1,[args])=LegalInterfaces then
print(`LegalInterfaces(m,n,LETTERr): Given positive integers`):
print(`m and n and a letter LETTERr with n boards, finds all the `):
print(`legal interfaces with the LETTERr on the right`):
print(`For example, try: LegalInterfaces(2,3,{{1},{2,3}}); `):
fi:
 
if nops([args])=1 and op(1,[args])=PreUmbra then
print(`PreUmbra(m,n,x,y,q,a,b,LETTERr,A): `):
print(`Given a letter from the right, LETTERr, with`):
print(`n boards, and an integer, m representing the max. number of boards`):
print(`of the letter on the left, and a symbol A, finds the`):
print(`pre-umbra (i.e. the image of a generic monomial`):
print(`x[1]^a[1]*y[1]^b[1]*...*y[n-1]^b[n-1]*x[m]^a[n]`):
print(`as a linear combination of A[LETTER], where LETTER`):
print(`denots m-board letters`):
print(`For example, try: PreUmbra(2,2,x,y,q,a,b,{{1},{2}},A); `):
fi:
 
if nops([args])=1 and op(1,[args])=PzA then
print(`PzA(z,A): given a list of variables, z, (of size r, say) and`):
print(`a symbolic integer A, finds the expression`):
print(`for the sum of z[1]^a1*...*z[r]^ar over all`):
print(`r-tuples of positive integers such that a1+...+ar=A`):
print(` For example, try PzA([x,y],B); `):
fi:
 
if nops([args])=1 and op(1,[args])=ToUmbra then
print(` ToUmbra(poly1,xlist,alist): given a polynomial, poly1, in the`):
print(` variables xlist with exponents that are affine-linear`):
print(` expressions in the discrete variables in the`):
print(` list alist, and in the variables of alist themselves`):
print(`outputs the corresponding Umbra, such that`):
print(` poly1 is the image, under that umbra, of the`):
print(` generic monomial y[1]^alist[1]*...*y[k]^alist[k],`):
print(`where k=nops(alist), and y[1], ..., y[k] are generic`):
print(` continuous variables that correspond to the disrete variables`):
print(` in alist`):
print(` The format of the output is a set, each element of which`):
print(` is a list of 3-elements`):
print(` [FRONT, Diffs,SubsList], where the FRONT `):
print(` is a rational function in whatever, Diffs is a list of `):
print(` integers of the same length as alist, and SubsList is `):
print(` the list of substitutions that the continuous variables `):
print(` that correspond to alist have to be substituted by `):
print(` For example ToUmbra(a*x^b+b*y^a,[x,y],[a,b]) should yield `):
print(`  {[1, [1, 0], [1, x]], [1, [0, 1], [y, 1]]} `):
fi:
 
 
 
end:
 
 
 
 
#NCP(kv): The set of Non-crossing partitions of a list kv
 
NCP:=proc(kv)
local kv1,la,gu,i,ru,kv2,mu,j,GU,nc:
with(combinat):
 
if nops(kv)=1 then
RETURN({{kv}}):
fi:
 
GU:={{kv}}:
 
la:=op(1,kv):
kv1:=kv minus {la} :
gu:=powerset(kv1) minus {kv1}:
for i from 1 to nops(gu) do
ru:=op(i,gu) union {la}:
kv2:=kv minus ru:
 
mu:=NCP(kv2):
 
for j from 1 to nops(mu) do
nc:=op(j,mu):
 
if not CRoss(nc,ru) then
GU:=GU union {nc union {ru}}:
fi:
 
od:
od:
GU:
 
end:
 
 
 
 
#CRoss(setP,kv): Checks whether a set partition is crossing with a set kv
CRoss:=proc(setP,kv)
local set1,i:
 
for i from 1 to nops(setP) do
 set1:=op(i,setP):
if Cross(set1,kv) then
 RETURN(true):
fi:
od:
false:
end:
 
 
 
#Cross(set1,set2): when two sets cross
Cross:=proc(set1,set2)
local a,b,c,d,i1,i2,j1,j2:
 
for i1 from 1 to nops(set1) do
a:=op(i1,set1):
for i2 from i1+1 to nops(set1) do
b:=op(i2,set1):
for j1 from 1 to nops(set2) do
c:=op(j1,set2):
for j2 from j1+1 to nops(set2) do
d:=op(j2,set2):
if cross({a,b},{c,d}) then
RETURN(true):
fi:
od:od:od:od:
false:
end:
 
 
#cross(set1,set2): when two 2-sets (say {a,b} and {c,d})
#cross i.e. a<c<b<d
cross:=proc(set1,set2)
local a,b,c,d:
a:=min(op(set1)):b:=max(op(set1)):
c:=min(op(set2)):d:=max(op(set2)):
 
if (a<c and c<b and b<d) or (c<a and a<d and d<b) then
RETURN(true):
else
RETURN(false):
fi:
end:
 
# Alphabet(L): The alphabet of 2D-lattice animals that have
# the property that each vertical cross-section has at most
# L boards (and hence L-1 gaps)
# 
Alphabet:=proc(L)
local gu,i,j:
option remember:
gu:={}:
for i from 1 to L do
gu:=gu union NCP({seq(j,j=1..i)}):
od:
gu:
end:
 
 
# Alphabet0(L): The alphabet of 2D-lattice animals that have
# the property that each vertical cross-section has 
# L boards (and hence L-1 gaps)
# 
Alphabet0:=proc(L)
local j:
NCP({seq(j,j=1..L)}):
end:
 
 
#IncSeq(n,K): All weakly increasing sequences of positive
#integers<=K of length n 
IncSeq:=proc(n,K)
local mu,gu,i,mu1,j:
if n=0 then
RETURN({[]}):
fi:
 
if n=1 then
RETURN({seq([i],i=1..K)}):
fi:
 
mu:=IncSeq(n-1,K):
gu:={}:
 
for i from 1 to nops(mu) do
mu1:=op(i,mu):
 
for j from mu1[n-1] to K do
gu:=gu union {[op(mu1),j]}:
od:
od:
gu:
end:
 
#Interfaces(m,n): All possible interfaces between an
#n-board animal-letter from the right and an m-board
#animal letter on its left
Interfaces:=proc(m,n)
local gu,mu,i,gu1,mu1,j:
 
mu:=IncSeq(2*m,2*n+1):
 
gu:={}:
 
for i from 1 to nops(mu) do
mu1:=op(i,mu):
gu1:=[seq([mu1[2*j-1],mu1[2*j]],j=1..m)]:
gu:=gu union {gu1}:
od:
 
gu:
 
end:
 
 
#DoesTouch(int1,interface1): given the number of the board
#(from bottom to top), int1, and an interface that defines
#a letter to the left, returns true if that interval
#gets touched by one of  the boards of the interface
DoesTouch:=proc(int1,interface1)
local i,mu:
for i from 1 to nops(interface1) do
mu:=op(i,interface1):
if mu[1]<=2*int1 and 2*int1<=mu[2] then
 RETURN(true):
fi:
od:
false:
end:
 
#DoesCompTouch(comp1,interface1): given a component of a letter
#comp1, and an interface that defines
#a letter to the left, returns true if that component
#gets touched by one of  the boards of the interface
DoesCompTouch:=proc(comp1,interface1)
local i:
 
for i from 1 to nops(comp1) do
 if DoesTouch(comp1[i],interface1) then
  RETURN(true):
 fi:
od:
 
false:
 
end:
 
 
#IsLegalInfc(LETTER1,interface1): given a letter LETTER1
#and an interface, interface1, that defines
#a letter to the left, returns true if the interface
#is a legal continuation
 
IsLegalInfc:=proc(LETTER1,interface1)
local i,comp1:
 
for i from 1 to nops(LETTER1) do
comp1:=op(i,LETTER1):
 
if not DoesCompTouch(comp1,interface1) then
 RETURN(false):
fi:
 
od:
 
true:
 
end:
 
#BoardTouchComp(board1,comp1): Given one board of an interface
#and one component of comp1, of a letter, returns true
#iff board1 touches the component comp1
BoardTouchComp:=proc(board1,comp1)
local i:
for i from 1 to nops(comp1) do
 
 if board1[1]<=2*comp1[i] and 2*comp1[i]<=board1[2] then
   RETURN(true):
 fi:
 
od:
 
false:
 
end:
 
#TouchingComps(board1,LETTER): Given a  letter LETTER,
#and a board, board1, finds the set of components of LETTER
#that touch board1
TouchingComps:=proc(board1,LETTER)
local i,gu:
gu:={}:
 
for i from 1 to nops(LETTER) do
 
if BoardTouchComp(board1,LETTER[i]) then
 gu:=gu union {LETTER[i]}:
fi:
od:
gu:
end:
 
########KAN
 
#ConComp(graph1,i): Given a vertex i in a graph1
#finds all the vertices in its connected component
ConComp:=proc(graph1,i)
local mu0,mu1,mu2,j:
mu0:={}:
mu1:={i}:
 
while mu0<>mu1 do
mu2:=mu1:
 
for j from 1 to nops(mu1) do
mu2:=mu2 union graph1[mu1[j]]:
od:
 
mu0:=mu1:
mu1:=mu2:
od:
mu2:
end:
 
#ConComps(graph1,V): Given a vertex set V and a graph1
#finds the connected components
ConComps:=proc(graph1,V)
local gu,i,V1:
V1:=V:
 
 
gu:={}:
 
while V1<>{} do
i:=V1[1]:
gu:=gu union {ConComp(graph1,i)}:
V1:=V1 minus ConComp(graph1,i):
od:
gu:
end:
 
 
 
#LeftLetter(Intfc1,LETTERr): Given an interface
#Intfc1, and a letter to its right, LETTERr,
#finds the letter than results on the left
#
LeftLetter:=proc(Intfc1,LETTERr)
local i,j,gu:
 
for i from 1 to nops(Intfc1) do
gu[i]:={}:
od:
 
for i from 1 to nops(Intfc1) do
for j from i+1 to nops(Intfc1) do
 
if TouchingComps(Intfc1[i],LETTERr) intersect
             TouchingComps(Intfc1[j],LETTERr)<>{}  then
gu[i]:=gu[i] union {j}:
gu[j]:=gu[j] union {i}:
fi:
od:
od:
ConComps(gu,{seq(i,i=1..nops(Intfc1))}):
end:
 
 
#LegalInterfaces(m,n,LETTERr): Given positive integers
#m and n and a letter LETTERr with n boards, finds all the 
#legal interfaces with the LETTERr on the right
LegalInterfaces:=proc(m,n,LETTERr)
local gu,mu,i:
 
mu:=Interfaces(m,n):
 
gu:={}:
 
for i from 1 to nops(mu) do
if  IsLegalInfc(LETTERr,mu[i]) then
gu:=gu union {mu[i]}:
fi:
od:
gu:
end:
 
#GapInterfaces0(Intfc): Need for GapInterfaces
#in terms of a list of intervals, finds the set of
#compatible interfaces for the gaps
GapInterfaces0:=proc(Intfc)
local Intfc1,i,a,b,n,j,gu,mu,lu:
 
n:=nops(Intfc):
if nops(Intfc)=1 then
a:=Intfc[1][1]:
if a=1 then
RETURN({[[1,1]]}):
else
RETURN({[[1,a]],[[1,a-1]]}):
fi:
fi:
 
Intfc1:=[op(1..n-1,Intfc)]:
 
mu:=GapInterfaces0(Intfc1):
 
a:=Intfc[n-1][2]:
b:=Intfc[n][1]:
 
if b=a then
lu:={[a,b]}:
elif b=a+1 then
lu:={[a,a],[a,b],[b,b]}:
elif b>a+1 then
lu:={[a,b],[a,b-1],[a+1,b],[a+1,b-1]}:
else
ERROR(`Bad input`):
fi:
 
gu:={}:
 
for i from 1 to nops(mu) do
for j from 1 to nops(lu) do
gu:=gu union {[op(mu[i]),lu[j]]}:
od:
od:
gu:
end:
 
#GapInterfaces(Intfc,n): Given an interface for the boards
#in terms of a list of intervals, finds the set of
#compatible interfaces for the gaps
#where n is the number of boards of the right letter
GapInterfaces:=proc(Intfc,n)
local gu,mu,lu,b,i,j:
mu:=GapInterfaces0(Intfc):
b:=Intfc[nops(Intfc)][2]:
 
if b=2*n+1 then
lu:={[b,b]}:
else
lu:={[b,2*n+1],[b+1,2*n+1]}:
fi:
 
 
gu:={}:
 
for i from 1 to nops(mu) do
for j from 1 to nops(lu) do
gu:=gu union {[op(mu[i]),lu[j]]}:
od:
od:
gu:
end:
 
 
#PzA(z,A): given a list of variables, z, (of size r, say) and
#a symbolic integer A, finds the expression
#for the sum of z[1]^a1*...*z[r]^ar over all
#r-tuples of positive integers such that a1+...+ar=A
PzA:=proc(z,A)
local i,r,z1,B,mu:
option remember:
r:=nops(z):
 
if r=1 then
RETURN(z[1]^A):
fi:
 
z1:=[op(2..r,z)]:
mu:=PzA(z1,B):
 
normal(simplify(normal(sum(z[1]^i*subs(B=A-i,mu),i=1..A-1)))):
end:
 
 
#APU(Intfc1,Gintfc1,m,n,x,y,q,a,b): Given an interface, Intfc1 for the
#m boards, and an interface Gintfc1, for the m+1 gaps,
#returns the Atomic Pre Umbra, corresponding to these
#Interfaces, i.e. the image of x[1]^a[1]*y[1]^b[1]*x[2]^a[2]...
#the right letter has n boards
APU:=proc(Intfc1,Gintfc1,m,n,x,y,q,a,b)
local A,j,i,gu,mu,gu1,gu2:
 
for i from 1 to 2*n+1 do
A[i]:={}:
od:
 
for i from 1 to 2*n+1 do
for j from 1 to m do
if Intfc1[j][1]<=i and i<=Intfc1[j][2] then
A[i]:=A[i] union {x[j]}:
fi:
od:
od:
 
 
for i from 1 to 2*n+1 do
for j from 1 to m+1 do
 
if Gintfc1[j][1]<=i and i<=Gintfc1[j][2] then
 A[i]:=A[i] union {y[j-1]}:
fi:
 
od:
od:
 
gu1:=1:
 
mu:=A[1] minus {y[0]}:
for i from 1 to nops(mu) do
gu1:=gu1*mu[i]/(1-mu[i]):
od:
 
 
gu2:=1:
 
mu:=A[2*n+1] minus {y[m]}:
for i from 1 to nops(mu) do
gu2:=gu2*mu[i]/(1-mu[i]):
od:
gu:=gu1*gu2:
 
for i from 1 to n do 
mu:=convert(A[2*i],list):
gu:=gu*PzA(mu,a[i]):
od:
 
for i from 1 to n-1 do 
mu:=convert(A[2*i+1],list):
gu:=gu*PzA(mu,b[i]):
od:
 
 
gu:=normal(subs({y[0]=1},gu)):
 
if normal(subs(y[m]=1,denom(gu)))<>0 then
gu:=normal(subs(y[m]=1,gu)):
else
gu:=
normal(subs(y[m]=1,normal(diff(numer(gu),y[m])/diff(denom(gu),y[m])))):
fi:
 
for i from 1 to m do
 gu:=subs(x[i]=q*x[i],gu):
od:
gu:=normal(simplify(expand(gu),symbolic)):
gu:
end:
 
 
#############Programs Burrowed from ROTA###############
# ToUmbra(poly1,xlist,alist): given a polynomial, poly1, in the
# variables xlist with exponents that are affine-linear
# expressions in the discrete variables in the
# list alist, and in the variables of alist themselves
# outputs the corresponding Umbra, such that
# poly1 is the image, under that umbra, of the
# generic monomial y[1]^alist[1]*...*y[k]^alist[k],
# (where k=nops(alist), and y[1], ..., y[k] are generic
# continuous variables that correspond to the disrete variables
# in alist
# The format of the output is a set, each element of whcih
# is a list of 3-elements
# [FRONT, Diffs,SubsList], where the FRONT
# is a rational function in whatever, Diffs is a list of
# integers of the same length as alist, and SubsList is
# the list of substitutions that the continuous variables
# that correspond to alist have to be substituted by
# For example ToUmbra(a*x^b+b*y^a,[x,y],[a,b]) should yield
#  {[1, [1, 0], [1, x]], [1, [0, 1], [y, 1]]}
ToUmbra:=proc(poly1,xlist,alist)
local gu,i,poly2:
poly2:=expand(poly1):
 
if not type(poly2,`+`) then
RETURN({UmbralTerm(poly2,xlist,alist)}):
fi:
 
gu:={}:
 
for i from 1 to nops(poly2) do
 gu:=gu union {UmbralTerm(op(i,poly2),xlist,alist)}:
od:
 
SimplifyUmbra(gu):
end:
 
 
# UmbralTerm(mono,xlist,alist): given a monomial in the
# variables xlist with exponents that are affine-linear
# expressions in the discrete variables in the
# list alist, and in the variables of alist themselves
# outputs the corresponding term in the Umbra
# The format of the output is a list of 3-elements
# [FRONT, Diffs,SubsList], where the FRONT
# is a rational function in whatever, Diffs is a list of
# integers of the same length as alist, and SubsList is
# the list of substitutions that the continuous variables
# that correspond to alist have to be substituted by
#
UmbralTerm:=proc(mono,xlist,alist)
local mono1,FRONT, Diffs, SubsList,i,d1,gu:
mono1:=expand(mono):
gu:=hafokh(mono1,xlist,alist):
FRONT:=gu[2]:
SubsList:=gu[1]:
 
mono1:=expand(FRONT):
Diffs:=[]:
 
for i from 1 to nops(alist) do
d1:=degree(mono1,alist[i]):
mono1:=expand(normal(expand(mono1/alist[i]^d1))):
Diffs:=[op(Diffs),d1]:
od:
 
[mono1,Diffs,SubsList]:
end:
 
 
 
#hafokh(mono,xlist,alist): given a monomial mono, in the form
#product(x[i]^a[j]) outputs the list [p[1],...p[k]] and the
#left-over monomial, lu, such that
#such that it equals lu*p[1]^a[1]*...*p[k]^a[k] times
hafokh:=proc(mono,xlist,alist)
local mono1,i,resh,mu:
mono1:=expand(mono):
resh:=[]:
for i from 1 to nops(alist) do
 mu:=grab(mono1,alist[i]):
 resh:=[op(resh),mu[1]]:
 mono1:=mu[2]:
od:
 
resh,mono1:
end:
 
 
 
#grab(mono,a): given a monomial, mono, in the variables 
#
# exponent (discrete) variable, a, returns
#the largest monomial, lu, such that lu^a is a factor of mono
grab:=proc(mono,a)
local mono1,mono2,i,lu,mu,mu1,mu2,khe,mu11,mu12:
mono1:=expand(mono):
mono2:=expand(mono1):
lu:=1:
 
if type(mono1,`*`)  then
 
for i from 1 to nops(mono1) do
mu:=op(i,mono1):
if type(mu,`^`) then
 mu1:=op(1,mu):
 mu2:=op(2,mu):
  if type(mu1,`^`) then
    mu11:=op(1,mu1):
     if type(mu11, `^`) then
      ERROR(`I give up`):
     fi:
    mu12:=op(2,mu1):
   mu1:=mu11:
   mu2:=mu2*mu12:
  fi:
 
 
 khe:=coeff(expand(mu2),a,1):
 lu:=lu*mu1^khe:
 mono2:=simplify(mono2/mu1^(a*khe)):
fi:
od:
 
elif type(mono1,`^`)  then
 mu1:=op(1,mono1):
 mu2:=op(2,mono1):
  if type(mu1,`^`) then
    mu11:=op(1,mu1):
     if type(mu11, `^`) then
      ERROR(`I give up`):
     fi:
    mu12:=op(2,mu1):
   mu1:=mu11:
   mu2:=mu2*mu12:
  fi:
 
 khe:=coeff(mu2,a,1):
 lu:=lu*mu1^khe:
 mono2:=simplify(mono2/mu1^(a*khe)):
fi:
lu,simplify(expand(mono2),symbolic):
end:
 
 
#SimplifyUmbra(Umb): Given an Umbra, collects like terms
#and returns a simplified version
SimplifyUmbra:=proc(Umb)
local Umb1,gu,kv,i,mu,evar,F:
Umb1:=ConvertToUmbra2(Umb):
 
gu:=0:
kv:={}:
 
for i from 1 to nops(Umb1) do
 gu:=gu+Umb1[i][1]*F(Umb1[i][2]):
 kv:=kv union {Umb1[i][2]}:
od:
 
gu:=expand(gu):
 
mu:={}:
 
for i from 1 to nops(kv) do
evar:=op(i,kv):
mu:=mu union {[factor(normal(coeff(gu,F(evar)))),evar[1],evar[2]]}:
od:
 
mu:
 
end:
 
#ConvertToUmbra2(Umb): Given an umbra Umb in the 3-list-format
#i.e. as a set of 3-lists of the form
#[FRONT,Orders_Of_Derivaties,Substitutions]
#converts this to a set of 2-lists of the form 
#[FRONT,[Orders_Of_Derivaties,Substitutions]]
ConvertToUmbra2:=proc(Umb)
local Umb1,i,mu:
Umb1:={}:
 
for i from 1 to nops(Umb) do
 mu:=op(i,Umb):
 Umb1:=Umb1 union {[mu[1],[mu[2],mu[3]]]}:
od:
 
Umb1:
end:
#############End Programs Burrowed from ROTA###############
 
 
#APUk(Intfc1,m,n,x,y,q,a,b): Compound atomic pre-umbra
#corresponding to Interface Intfcl
APUk:=proc(Intfc1,m,n,x,y,q,a,b)
local mu1, gu1, j:
gu1:=0:
 
mu1:=GapInterfaces(Intfc1,n):
for j from 1 to nops(mu1) do
gu1:=gu1+APU(Intfc1,mu1[j],m,n,x,y,q,a,b):
od:
normal(simplify(gu1)):
end:
 
 
 
#PreUmbra(m,n,x,y,q,a,b,LETTERr,A): 
#Given a letter from the right, with
#n boards, and an integer, m representing the maximal
#number of boards
#of the letter on the left, and a symbol A, finds the
#pre-umbra (i.e. the image of a generic monomial
#x[1]^a[1]*y[1]^b[1]*...*y[n-1]^b[n-1]*x[m]^a[n]
#as a linear combination of A[LETTER], where LETTER
#denots m-board letters
 
PreUmbra:=proc(m,n,x,y,q,a,b,LETTERr,A)
local gu,mu,i,LeftL,gu1,m1: 
gu:=0:
 
for m1 from 1 to m do
mu:=LegalInterfaces(m1,n,LETTERr):
 
for i from 1 to nops(mu) do
LeftL:=LeftLetter(mu[i],LETTERr):
gu1:=APUk(mu[i],m1,n,x,y,q,a,b):
gu:=gu+gu1*A[LeftL]:
od:
od:
gu:=expand(gu):
 
mu:=Alphabet(m):
gu1:=0:
for i from 1 to nops(mu) do
gu1:=gu1+normal(coeff(gu,A[mu[i]],1))*A[mu[i]]:
od:
 
gu1:
 
end:
 
 
 
#AniUmSc(m,x,y,q): An Umbral Scheme, i.e. a list 
#with four entries:
#a list of lists of Umbra, followed by the list of
#arguments of the corresponding functions
#for m-board lattice animals
#followed by the initial values of the corresponding function
#followed by the set of final letters (or states)
AniUmSc:=proc(m,x,y,q)
local mu,gu,n,a,b,A,lu,T,S,mu1,j,i1,xlist,alist,i,tu,gu1,U,zu,yu,zu1:
mu:=Alphabet(m):
zu:={}:
gu:=[]:
for n from 1 to m do
mu1:=Alphabet0(n):
 
for i from 1 to nops(mu1) do
xlist:=[x[1],seq(op([y[i1],x[i1+1]]),i1=1..n-1)]:
alist:=[a[1],seq(op([b[i1],a[i1+1]]),i1=1..n-1)]:
S[mu1[i]]:=xlist:
 
if nops(mu1[i])=n then
 U[mu1[i]]:=Ini(x,y,q,n):
else
 U[mu1[i]]:=0:
fi:
 
if nops(mu1[i])=1 then
zu:=zu union {mu1[i]}:
fi:
 
lu:=PreUmbra(m,n,x,y,q,a,b,mu1[i],A):
 
for j from 1 to nops(mu) do
T[mu1[i],mu[j]]:=ToUmbra(coeff(lu,A[mu[j]],1),xlist,alist):
od:
od:
od:
 
gu:=[]:
tu:=[]:
 
mu:=convert(mu,list):
 
yu:=[]:
zu1:={}:
for i from 1 to nops(mu) do
 tu:=[op(tu),S[mu[i]]]:
 yu:=[op(yu),U[mu[i]]]:
 if member(mu[i],zu) then
  zu1:=zu1 union {i}:
 fi:
 
 
 gu1:=[]:
 
 for j from 1 to nops(mu) do
  gu1:=[op(gu1),T[mu[j],mu[i]]]:
 od:
 
 gu:=[op(gu),gu1]:
od:
 
[zu1,gu,yu,tu]:
 
end:
 
#Ini(x,y,q,m): The initial weight of a leftmost letter
#with m boards according to x (for the boards), y (for the gaps)
#and q (for the length of the boards)
Ini:=proc(x,y,q,m)
local gu,i:
gu:=1:
 
for i from 1 to m do
 gu:=gu*q*x[i]/(1-q*x[i]):
od:
 
for i from 1 to m-1 do
 gu:=gu*y[i]/(1-y[i]):
od:
gu:
end:
 
 
##############The Pre-Computed Animal Schemes
 
AUS1:=
[{1}, [[{[1/(-1+q*x[1])**2*q**2*x[1]**2, [0], [1]], [-q*x[1]/(-1+q*x[1]), [1]
, [1]]}]], [q*x[1]/(1-q*x[1])], [[x[1]]]]:
AUS2:=
[{1, 3}, [[{[q**3*x[2]**2*x[1]/(-y[1]+q*x[2])/(-1+y[1])/(-1+q*x[2])/(-1+q*x[1])
, [0, 0, 0], [y[1], q*x[2], 1]], [q**2*x[2]*x[1]/(-1+y[1])**2/(-1+q*x[2])/(-1+q
*x[1]), [0, 0, 0], [y[1], y[1], y[1]]], [q**3*x[2]*x[1]**2/(-1+y[1])/(q*x[1]-y[
1])/(-1+q*x[2])/(-1+q*x[1]), [0, 0, 0], [1, q*x[1], y[1]]], [q**2*x[2]*x[1]*y[1
]**2/(-y[1]+q*x[2])/(-1+y[1])**2/(q*x[1]-y[1]), [0, 0, 0], [1, y[1], 1]], [-q**
2*x[2]*y[1]*x[1]/(-y[1]+q*x[2])/(-1+y[1])**2/(-1+q*x[1]), [0, 0, 0], [y[1], y[1
], 1]], [-q**2*x[2]*y[1]*x[1]/(-1+y[1])/(-1+q*x[2])/(-1+q*x[1]), [1, 0, 0], [1
, 1, 1]], [-q**2*x[2]*y[1]*x[1]/(-1+y[1])/(-1+q*x[2])/(-1+q*x[1]), [0, 0, 1], [
1, 1, 1]], [q**2*x[2]*x[1]/(-1+y[1])**2/(-1+q*x[2])/(-1+q*x[1]), [0, 0, 0], [1
, 1, y[1]]], [2*q**2*x[2]*y[1]*x[1]*(y[1]-2)/(-1+y[1])**2/(-1+q*x[2])/(-1+q*x[1
]), [0, 0, 0], [1, 1, 1]], [-q**2*x[2]*y[1]*x[1]/(-1+y[1])**2/(q*x[1]-y[1])/(-1
+q*x[2]), [0, 0, 0], [1, y[1], y[1]]], [-x[1]**2*y[1]*x[2]*q**2*(q*x[1]-q*x[2]+
1-y[1])/(-1+y[1])/(q*x[1]-y[1])/(x[1]-x[2])/(-1+q*x[2])/(-1+q*x[1]), [0, 0, 0]
, [1, q*x[1], 1]], [q**2*x[2]*x[1]/(-1+y[1])**2/(-1+q*x[2])/(-1+q*x[1]), [0, 0
, 0], [y[1], 1, 1]], [-x[1]*y[1]*x[2]**2*q**2*(q*x[1]-q*x[2]-1+y[1])/(-y[1]+q*x
[2])/(-1+y[1])/(x[1]-x[2])/(-1+q*x[2])/(-1+q*x[1]), [0, 0, 0], [1, q*x[2], 1]]}
, {[q**2*x[2]*x[1]*y[1]/(-1+q*x[1])/(-1+q*x[2])**2/(-y[1]+q*x[2]), [0, 0, 0], [
q*x[2], q*x[2], 1]], [-q**3*x[2]*x[1]**2/(-1+y[1])/(q*x[1]-y[1])/(-1+q*x[2])/(-
1+q*x[1]), [0, 0, 0], [1, q*x[1], y[1]]], [q**3*x[2]**2*x[1]*y[1]/(-1+q*x[1])/(
-1+q*x[2])**2/(-1+y[1]), [0, 0, 0], [1, q*x[2], 1]], [-q**3*x[2]**2*x[1]/(-y[1]
+q*x[2])/(-1+y[1])/(-1+q*x[2])/(-1+q*x[1]), [0, 0, 0], [y[1], q*x[2], 1]], [q**
3*x[2]*x[1]**2*y[1]/(-1+q*x[1])**2/(-1+q*x[2])/(-1+y[1]), [0, 0, 0], [1, q*x[1]
, 1]], [q**2*x[2]*x[1]*y[1]/(-1+q*x[1])**2/(-1+q*x[2])/(q*x[1]-y[1]), [0, 0, 0]
, [1, q*x[1], q*x[1]]]}, {[-q**2*x[2]*y[1]*x[1]/(-1+y[1])/(-1+q*x[2])/(-1+q*x[1
]), [1], [1]], [q**2*x[2]*y[1]*x[1]*(y[1]-2)/(-1+y[1])**2/(-1+q*x[2])/(-1+q*x[1
]), [0], [1]], [q**2*x[2]*x[1]/(-1+y[1])**2/(-1+q*x[2])/(-1+q*x[1]), [0], [y[1]
]]}], [{[-2*q**3*x[2]**2*x[1]/(-y[1]+q*x[2])/(-1+y[1])/(-1+q*x[2])/(-1+q*x[1])
, [0, 0, 0], [y[1], q*x[2], 1]], [2*q**2*x[2]*y[1]*x[1]/(-1+y[1])**2/(q*x[1]-y[
1])/(-1+q*x[2]), [0, 0, 0], [1, y[1], y[1]]], [-2*q**3*x[2]*x[1]**2/(-1+y[1])/(
q*x[1]-y[1])/(-1+q*x[2])/(-1+q*x[1]), [0, 0, 0], [1, q*x[1], y[1]]], [-2*q**2*x
[2]*x[1]*y[1]**2/(-y[1]+q*x[2])/(-1+y[1])**2/(q*x[1]-y[1]), [0, 0, 0], [1, y[1]
, 1]], [2*q**2*x[2]*y[1]*x[1]*(-3*q*x[2]+q*x[2]*y[1]+y[1]*q*x[1]-3*q*x[1]-2*y[1
]+4+2*q**2*x[2]*x[1])/(-1+q*x[1])**2/(-1+q*x[2])**2/(-1+y[1])**2, [0, 0, 0], [1
, 1, 1]], [2*q**2*x[2]*y[1]*x[1]/(-y[1]+q*x[2])/(-1+y[1])**2/(-1+q*x[1]), [0, 0
, 0], [y[1], y[1], 1]], [(-q**2*x[2]**2+q**2*x[2]*x[1]+y[1]*q*x[1]-2*q*x[1]+q*x
[2]*y[1]+2-2*y[1])*q**2*x[2]**2*y[1]*x[1]/(x[1]-x[2])/(-y[1]+q*x[2])/(-1+q*x[2]
)**2/(-1+q*x[1])/(-1+y[1]), [0, 0, 0], [1, q*x[2], 1]], [(q**2*x[1]**2-q**2*x[2
]*x[1]-q*x[2]*y[1]+2*q*x[2]-y[1]*q*x[1]-2+2*y[1])*q**2*x[2]*y[1]*x[1]**2/(x[1]-
x[2])/(-1+q*x[2])/(q*x[1]-y[1])/(-1+q*x[1])**2/(-1+y[1]), [0, 0, 0], [1, q*x[1]
, 1]], [-2*q**2*x[2]*x[1]/(-1+y[1])**2/(-1+q*x[2])/(-1+q*x[1]), [0, 0, 0], [1, 
1, y[1]]], [-2*q**2*x[2]*x[1]/(-1+y[1])**2/(-1+q*x[2])/(-1+q*x[1]), [0, 0, 0], 
[y[1], 1, 1]], [-2*q**2*x[2]*x[1]/(-1+y[1])**2/(-1+q*x[2])/(-1+q*x[1]), [0, 0, 
0], [y[1], y[1], y[1]]]}, {[q**2*x[2]*x[1]/(-1+y[1])**2/(-1+q*x[2])/(-1+q*x[1])
, [0, 0, 0], [y[1], y[1], y[1]]], [q**2*x[2]*x[1]*y[1]**2/(-y[1]+q*x[2])/(-1+y[
1])**2/(q*x[1]-y[1]), [0, 0, 0], [1, y[1], 1]], [-q**2*x[2]*y[1]*x[1]/(-y[1]+q*
x[2])/(-1+y[1])**2/(-1+q*x[1]), [0, 0, 0], [y[1], y[1], 1]], [-q**2*x[2]*y[1]*x
[1]/(-1+y[1])**2/(q*x[1]-y[1])/(-1+q*x[2]), [0, 0, 0], [1, y[1], y[1]]], [-x[1]
**2*y[1]*x[2]*q**2*(q*x[1]-q*x[2]+1-y[1])/(-1+y[1])/(q*x[1]-y[1])/(x[1]-x[2])/(
-1+q*x[2])/(-1+q*x[1]), [0, 0, 0], [1, q*x[1], 1]], [-x[1]*y[1]*x[2]**2*q**2*(q
*x[1]-q*x[2]-1+y[1])/(-y[1]+q*x[2])/(-1+y[1])/(x[1]-x[2])/(-1+q*x[2])/(-1+q*x[1
]), [0, 0, 0], [1, q*x[2], 1]], [2*q**3*x[2]**2*x[1]/(-y[1]+q*x[2])/(-1+y[1])/(
-1+q*x[2])/(-1+q*x[1]), [0, 0, 0], [y[1], q*x[2], 1]], [-q**2*x[2]*x[1]*y[1]/(-
1+q*x[1])**2/(-1+q*x[2])/(q*x[1]-y[1]), [0, 0, 0], [1, q*x[1], q*x[1]]], [2*q**
3*x[2]*x[1]**2/(-1+y[1])/(q*x[1]-y[1])/(-1+q*x[2])/(-1+q*x[1]), [0, 0, 0], [1, 
q*x[1], y[1]]], [-q**2*x[2]*x[1]*y[1]/(-1+q*x[1])/(-1+q*x[2])**2/(-y[1]+q*x[2])
, [0, 0, 0], [q*x[2], q*x[2], 1]]}, {[q**2*x[2]*y[1]*x[1]*(-3*q*x[2]+q*x[2]*y[1
]+y[1]*q*x[1]-3*q*x[1]-2*y[1]+4+2*q**2*x[2]*x[1])/(-1+q*x[1])**2/(-1+q*x[2])**2
/(-1+y[1])**2, [0], [1]], [-2*q**2*x[2]*x[1]/(-1+y[1])**2/(-1+q*x[2])/(-1+q*x[1
]), [0], [y[1]]]}], [{[2/(-1+q*x[1])**2*q**2*x[1]**2, [0, 0, 0], [1, 1, 1]], [-
1/(-1+q*x[1])**2*q**2*x[1]**2, [0, 0, 0], [1, q*x[1], 1]], [-q*x[1]/(-1+q*x[1])
, [0, 0, 1], [1, 1, 1]], [-q*x[1]/(-1+q*x[1]), [1, 0, 0], [1, 1, 1]]}, {[1/(-1+
q*x[1])**2*q**2*x[1]**2, [0, 0, 0], [1, q*x[1], 1]]}, {[1/(-1+q*x[1])**2*q**2*x
[1]**2, [0], [1]], [-q*x[1]/(-1+q*x[1]), [1], [1]]}]], [0, q**2*x[1]/(1-q*x[1])
*x[2]/(1-q*x[2])*y[1]/(-y[1]+1), q*x[1]/(1-q*x[1])], [[x[1], y[1], x[2]], [x[1]
, y[1], x[2]], [x[1]]]]:
 
 
##############End Pre-Computed Animal Scheme###############
 
 
 
 
 
# ApplyUmbralMatrix(UmbraM,Fvec,ylistvec): Given an Umbral matrix
# UmbraM (in terms of a list of lists), a vector of expressions,
# Fvec, and a vector to variable-lists corresponding to the
# components of UmbraM and Fvec
#
ApplyUmbralMatrix:=proc(UmbraM,Fvec,ylistvec)
local gu,mu,i,j,k:
 
k:=nops(UmbraM):
 
if k<>nops(Fvec) or k<>nops(ylistvec) then
ERROR(`Bad input`):
fi:
 
gu:=[]:
 
for i from 1 to k do
mu:=0:
for j from 1 to k do
mu:=normal(mu+ApplyUmbra(UmbraM[i][j],Fvec[j],ylistvec[j])):
od:
gu:=[op(gu),mu]:
od:
gu:
end:
 
 
# ApplyUmbra(Umbra1,f,ylist): given an umbra, Umbra1, applies
# it to the expression f in the variables in ylist
ApplyUmbra:=proc(Umbra1,f,ylist)
local i,gu:
gu:=0:
 
for i from 1 to nops(Umbra1) do
 gu:=normal(gu+ApplyUmbraTerm(Umbra1[i],f,ylist)):
od:
gu:
end:
 
 
# ApplyUmbraTerm(Umbra1,f,ylist): given an umbral term, Umbra1, applies
# it to the expression f in the variables in ylist
ApplyUmbraTerm:=proc(Umbra1,f,ylist)
local i, FRONT,Diffs,SubsList,gu,bu:
 
FRONT:=Umbra1[1]:
Diffs:=Umbra1[2]:
SubsList:=Umbra1[3]:
 
if nops(Diffs)<>nops(ylist) then
 ERROR(`Diffs and ylist should have the same length`):
fi:
 
 
if nops(SubsList)<>nops(ylist) then
 ERROR(`SubsList and ylist should have the same length`):
fi:
 
gu:=f:
 
for i from 1 to nops(ylist) do
gu:=normal(xdiff1(gu,ylist[i],Diffs[i])):
od:
 
bu:={}:
for i from 1 to nops(ylist) do
bu:=bu union {ylist[i]=SubsList[i]}:
od:
gu:=subs(bu,gu):
 
normal(FRONT*gu):
 
end:
 
 
# xdiff1(f,x,i): given an expression f, and a variable x,
# and an integer i, finds (xD_x)^i f
# 
xdiff1:=proc(f,x,i)
local gu,j:
 
if i=0 then
 RETURN(f):
fi:
 
gu:=f:
 
for j from 1 to i do
gu:=x*diff(gu,x):
od:
 
gu:
 
end:
 
 
 
 
# ApplyUmSc(UmSch,q,n,var): Given an Umbral Scheme, UmSch, finds
# the generating functions of creatures up to the wt. q^n
#followed by substition of all x and y variables to 1 and summing
#over leftmost letters
ApplyUmSc:=proc(UmSch,q,n,var)
local UM,INI,ylistvec,lu,gu,i,lu1,i1,kv,ku,fu,Sid,j,mu,gc:
if not type(n,integer) or not n>0 then
ERROR(`n>0`):
fi:
 
 
kv:=UmSch[1]:
UM:=UmSch[2]:
INI:=UmSch[3]:
ylistvec:=UmSch[4]:
 
lu:=[]:
 
for i from 1 to nops(INI) do
lu:=[op(lu),SerToPol(INI[i],q,0)]:
od:
 
if n=0 then
RETURN(lu):
fi:
 
Sid:=[]:
 
for i from 1 to n do
 
 lu1:=ApplyUmbralMatrix(UM,lu,ylistvec):
 gu:=[]:
 
for i1 from 1 to nops(lu1) do
 gu:=[op(gu),expand(SerToPol(INI[i1],q,i)+SerToPol(lu1[i1],q,i))]:
od:
lu:=gu:
 
fu:=[]:
for i1 from 1 to nops(lu) do
fu:=[op(fu),coeff(expand(lu[i1]),q,i)]:
od:
 
mu:=0:
 
for j from 1 to nops(kv) do
 mu:=mu+fu[kv[j]]:
od:
 
 
ku:={seq(var[gc]=1,gc=1..nops(var))}:
mu:=subs(ku,mu):
 
Sid:=[op(Sid),factor(normal(mu))]:
 
lprint(Sid):
 
od:
 
RETURN(lu,Sid):
 
end:
 
 
#SerToPol(sidra,q,n): given a series, sidra, finds the polynomial
#consisting of the terms up to degree n
SerToPol:=proc(sidra,q,n)
local lu,i,gu:
lu:=taylor(sidra,q=0,n+1):
gu:=0:
for i from 0 to n do
gu:=gu+coeff(lu,q,i)*q^i:
od:
gu:
end: