######################################################################
##MAYLIS: Save this file as MAYLIS. To use it, stay in the                 #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read MAYLIS : <Enter>                                 #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Temple University ,                   #
#zeilberg@math.temple.edu.                                           # 
#######################################################################

#Created: Jan.-March 2001
#This version: March 5, 2001
#MAYLIS: A Maple package that is targeted to treat classes of 
#self-avoiding polygons with at most 2 horizontal edges per
#vertical cross-section, 
#using the Umbral Transfer Matrix Theorem
#It is one of the packages accompanies Doron Zeilberger's article:
#The Umbral Transfer Matrix Method. IV. Counting Self-Avoiding
#Polygons and Walks
#Please report bugs to zeilberg@math.temple.edu



print(`MAYLIS: A Maple package that is targeted to treat classes of `):
print(`self-avoiding polygons with at most 2 horizontal edges per`):
print(`vertical cross-section, `):
print(`Equivalently, row-convex polyominoes according to perimeter`):
print(`using the Umbral Transfer Matrix Theorem`):
print(`It is one of the packages that accompany Doron Zeilberger's article:`):
print(`The Umbral Transfer Matrix Method. IV. Counting Self-Avoiding`):
print(`Polygons and Walks`):
lprint(``):
print(`Created: March 8, 2001`):
print(`This version: March 8, 2001`):
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 IMPORTANT procedures type ezra(), for help with`):
 print(`a specific procedure, type ezra(procedure_name)`):
 print(`For a list of ALL procedures that are specific to MAYLIS`):
print(`i.e. are not part of USAP, type ezra1();, for help with`):
 print(`a specific procedure, type ezra(procedure_name)`):

print(`SAPUS1 and SAPUS2 are the pre-computed Umbral Schemes for`):
print(`1-LR and 2-LR self-avoiding polygons`):
 print(``):


ezra1:=proc()
print(``):
if nops([args])=0  then
print(`contains procedures: Alphabet, ApplyUmSc, ApplyUmSc ` ):
print(` APU, Delest, Delest1, DoThemAll,DoThemAllT, DV, FindL, FindR, `):
print(` Followers ,FollowersR, FunEq, JoinLL, JoinRL, JoinRR , Maylis,`):
print(` OneStepAll, OutLetter, PreFollowers, `):
print(` PrePreFollowers, PreUmbra, ProveDelestViennot`):
print(` PzA , SAPUmSc , SAPUmScg, `):
print(` SAPUmScgR, SAP1UmScR, SSUE, SSUEx, StrictDV, ToUmbra `):
print( ` Xavier, Xavierx, YafeUm, YafeUmSc `):
fi:
end:
ezra:=proc()
print(``):
if nops([args])=0  then
print(`The main  procedures are`):
print(` Delest, DoThemAll,DoThemAllT, DV, FunEq, Maylis `):
print(` ProveDelestViennot, SAPUmScgR, SAP1UmScR, SSUE, SSUEx `):
print( ` Xavier, Xavierx`):
fi:

if nops([args])=1 and op(1,[args])=Alphabet then
print(` Alphabet(n,L,R): The alphabet of 2D-self-avoiding polygons that have`):
print(` the property that each vertical cross-section has at most`):
print(` n Left and n right parantheses `):
print(`For example, to get the alphabet for <=2 pairs do: Alphabet(2); `)
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])=ApplyUmScg then
print(` ApplyUmScg(UmSch,q,n,var): Given a Gen. 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(InpLet,OutPreLet,L,R,C,x,q,A):`):
print(` Given an input letter from the left`):
print(`InpLet, and an output Pre-letter from the right,`):
print(`OutPreLet, finds the Atomic Pre-Umbra`):
print(`Here x is the indexed variable, and q is the "length variable"`):
fi:

if nops([args])=1 and op(1,[args])=Delest then
print(`Delest(RsetL,q): The Generating Function for column-convex`):
print(`polyominoes (equivalenty, SAPs with at most two horizontal`):
print(`edges per vertical cross-section) according to parameter,`):
print(`where it firsts goes according to the`):
print(`interfaces in RsetL[1], then it either quits or starts`):
print(`going according to the`):
print(`interfaces in RsetL[2], then it quits, or starts going`):
print(`according to the interfaces in RsetL[3], and so on`):
fi:

if nops([args])=1 and op(1,[args])=Delest1 then
print(`Delest1(RsetL,q): The Generating Function for column-convex`):
print(`polyominoes (equivalenty, SAPs with at most two horizontal`):
print(`edges per vertical cross-section) according to parameter,`):
print(`where it firsts goes according to the`):
print(`interfaces in RsetL[1], then it starts`):
print(`going according to the`):
print(`interfaces in RsetL[2], then  starts going`):
print(`according to the interfaces in RsetL[3], and so on`):
print(`Unlike Delest, it is not allowed to quit in the middle`):
fi:

if nops([args])=1 and op(1,[args])=DoThemAll then
print(`DoThemAll(q): Tries to Find the generating funtion of`):
print(`All Consistent classes of Column-Convex polyominoes`):
print(`with the further restriction that the interfaces`):
print(`obey the set Rset, for all possible (inequivalent)`):
print(`subsets of the full gamut (for column-convex polyominoes)`):
fi:

if nops([args])=1 and op(1,[args])=DoThemAllT then
print(`DoThemAllT(TOT,q): Tries to Find the generating funtion of`):
print(`All Consistent classes of Column-Convex polyominoes`):
print(`with the further restriction that the interfaces`):
print(`obey the set Rset, for all possible (inequivalent)`):
print(`subsets of the set of interfaces TOT `):
fi:

if nops([args])=1 and op(1,[args])=DV then
print(`DV(q): The Delest-Viennot Generating Function for the`):
print(`number of convex polygons according to perimeter`):
fi:

if nops([args])=1 and op(1,[args])=FindL then
print(` FindL(w,j,L,R): Given a legal word w (in {L,R}), and`):
print(` a place j that is an R, find the location of its L-mate`):
fi:

if nops([args])=1 and op(1,[args])=FindR then
print(` FindR(w,i,L,R): Given a legal word w (in {L,R}), and`):
print(` a place i that is an L, find the location of its R-mate`):
fi:

if nops([args])=1 and op(1,[args])=Followers then
print(`Followers(LET,L,R,C,n): All the ways to follow`):
print(`the letter LET if the maximal number of pairs  of`):
print(`parantheses is n`):
fi:

if nops([args])=1 and op(1,[args])=FollowersR then
print(`FollowersR(LET,L,R,C,n,Rset): All the ways to follow`):
print(`the letter LET if the maximal number of pairs  of`):
print(`parantheses is n, where the extensions should be`):
print(`in the set Rset`):
fi:

if nops([args])=1 and op(1,[args])=FunEq then
print(`FunEq(SUE,A): Given a Simplified Umbral Scheme`):
print(`(with one function), and a symbol A, finds`):
print(`the Functional Equation in a Human Readable Form`):
fi:

if nops([args])=1 and op(1,[args])=JoinLL then
print(`JoinLL(LET1,ClosedI,i,L,R,C): Given a letter LET1,`):
print(` (phrased in terms of L and R and C(for closed) ), `):
print(` A set of intervals of the form [a,a+1] indicating`):
print(` "closed blocks", and a location i, for which`):
print(` such that LET1[1][i]=LET1[1][i+1]=L does the`):
print(` operation LLpRqR->CCpLqR and appends to ClosedI `):
print(` the interval [i,i+1]`):
fi:

if nops([args])=1 and op(1,[args])=JoinRL then
print(`JoinRL(LET1,ClosedI,j,L,R,C): Given a letter LET1,`):
print(` (phrased in terms of L and R), `):
print(` such that LET1[1][j-1]=R, LET1[1][j]=L does the`):
print(` operation pRLq->pCCq and appends [j,j+1] to `):
print(` on ClosedI `):
fi:

if nops([args])=1 and op(1,[args])=JoinRR then
print(`JoinRR(LET1,ClosedI,j,L,R,C): Given a letter LET1,`):
print(` (phrased in terms of L and R and C(for closed) ), `):
print(` A set of intervals of the form [a,a+1] indicating`):
print(` "closed blocks", and a location i, for which`):
print(` such that LET1[1][i]=LET1[1][i+1]=L does the`):
print(` operation LpLqRR->LpRqCC `):
print(` operation LLpRqR->CCpLqR and appends to ClosedI `):
print(` the interval [i,i+1] `):
fi:


if nops([args])=1 and op(1,[args])=Maylis then
print(`Maylis(x,q,Rset,A): `):
print(`The Functional Equation, in Human Form, for the`):
print(`generating function, let's call it A(x) for SAPS`):
print(`with exactly two horizontal edges in every vertical`):
print(`cross-section, with the interfaces restriced by the`):
print(`set of pairs Rset`):
print(`according to the weight q^parameter x^(width of last vertical wall)`):
fi:

if nops([args])=1 and op(1,[args])=OneStepAll then
print(`OneStepAll(LET1,ClosedI,L,R,C): Set of all pairs `):
print(` [LET1, ClosedI]`):
print(` that can be obtained `):
print(` from the letter LET1, with its accompanying ClosedI `):
print(` L an R denote ( and ) resp. and C is neutral `):
fi:

if nops([args])=1 and op(1,[args])=OutLetter then
print(`OutLetter(PreLet,L,R,C): Given a pre-letter`):
print(`[Word,ClosedI,FS,NewLR_assignment], finds`):
print(`the resulting letter that emerges on the right`):
fi:

if nops([args])=1 and op(1,[args])=PreFollowers then
print(`PreFollowers(LET,L,R,C): Given an Umbral letter in`):
print(`the SAP alphabet (a word in the alphabet {L,R})`):
print(`finds all the pre-followers`):
fi:

if nops([args])=1 and op(1,[args])=PrePreFollowers then
print(`PrePreFollowers(LET,L,R,C): Given an Umbral  letter LET`):
print(`and symbols L,R,C, finds the set of all pre-pre-followers`):
print(` which are pairs of the form [LET1,ClosedI]`):
fi:

if nops([args])=1 and op(1,[args])=PreUmbra then
print(`PreUmbra(LETTER,L,R,n,x,q,A,Z,HATKH,SOF)`):
print(`Given a letter LETTER, in terms of legal`):
print(`parantheses using L and R`):
print(`and an integer n specifying the max. number of`):
print(`LR-pairs, and a letter Z`):
print(`and a symbol A, finds the`):
print(`pre-umbra (i.e. the image of a generic monomial`):
print(`x[1]^A[1]*x[2]^A[2]*...*x[r]^A[r] where r is nops(LETTER)-1`):
print(`as a linear combination of Z[LETTER1], where LETTER1`):
print(`ranges over all followers of LETTER`):
print(` HATKH and SOF denote the INITIAL and TERMINAL letters, respectively`):
fi:

if nops([args])=1 and op(1,[args])=ProveDelestViennot then
print(`ProveDelestViennot(): Proves the Delest-Viennot Generating Function`):
print(`(equivalent to the explicit formula `):
print(`a(2n+8)=(2n+11)4^n-4(2n+1)!/n!^2,`):
print(`for n>=0, and a(4)=1, a(6)=2)`):
print(`It returns true if the formula is correct, and false otherwise`):
print(`Note that this furnishes a completely automatic RIGOROUS proof`):
print(`of this celebrated result`):
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`):
fi:

if nops([args])=1 and op(1,[args])=SAPUmSc then
print(`SAPUmSc(n,x,q): An Umbral Scheme, i.e. a list `):
print(`with four entries:`):
print(`a list of lists of Umbra, followed by the list of`):
print(`arguments of the corresponding functions`):
print(`followed by the initial values of the corresponding function`):
print(`followed by the set of final letters (or states)`):
print(`for self-avoiding polygons with the property`):
print(`that each vertical cross section has at most n LR-pairs`):
fi:

if nops([args])=1 and op(1,[args])=SAPUmScg then
print(`SAPUmScg(n,x,q,t,w): An Umbral Scheme, i.e. a list `):
print(`with four entries:`):
print(`a list of lists of Umbra, followed by the list of`):
print(`arguments of the corresponding functions`):
print(`followed by the initial values of the corresponding function`):
print(`followed by the set of final letters (or states)`):
print(`for self-avoiding polygons with the property`):
print(`that each vertical cross section has at most n LR-pairs`):
print(`here the weight is q^(perimeter)*t^(area)*w^(width)`):
fi:

if nops([args])=1 and op(1,[args])=SAPUmScgR then
print(`SAPUmScgR(n,x,q,t,w,Rset): `):
print(`Like SAPUmScgR(n,x,q,t,w), but only for SAPS for which the`):
print(`continuations belong to the set Rset`):
fi:

if nops([args])=1 and op(1,[args])=SAP1UmScR then
print(`SAP1UmScR(x,q,Rset): `):
print(`A simplified Umbral Scheme in the form`):
print(`[INI,Rota_Operator]`):
print(`for SAPS with only two horizontal edges per vertical `):
print(`cross-section, in which the follow-up rules`):
print(`are according to Rset, where 1 means up, 0 means level,`):
print(`and -1 means down`):
print(`For example SAP1UmScR(x,q,{[0,0],[0,1],[1,0],[1,1]}) `):
print(`gives the Umbral Scheme for suchs SAPS in which`):
print(`you never go down from a horizontal edge, but have`):
print(`a weakly upwards trend`):
fi:

if nops([args])=1 and op(1,[args])=SSUE then
print(`SSUE(SUE,q): Given a Simplified Umbral Equation, SUE, and`):
print(` a variable q, that is a linear combinatrion of`):
print(`A(x),A(q), and qA'(q), solves for A(q)`):
fi:

if nops([args])=1 and op(1,[args])=SSUEx then
print(`SSUEx(SUE,q,x): Given a Simplified Umbral Equation, SUE, and`):
print(` a variable q, that is a linear combinatrion of`):
print(`A(x),A(q), and qA'(q), solves for A(q)`):
print(`it also returns A(x)`):
fi:

if nops([args])=1 and op(1,[args])=StrictDV then
print(`DV(q): The Analog of the Delest-Viennot Generating Function for the`):
print(`number of Strictly convex polygons according to perimeter`):
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:

if nops([args])=1 and op(1,[args])=Xavier then
print(`Xavier(Rset,q): The Generating Function for column-convex`):
print(`polyominoes (equivalenty, SAPs with at most two horizontal`):
print(`edges per vertical cross-section) according to perimeter,`):
print(`where interfaces are constrained to belong to the set`):
print(`Rset, for example, try: Xavier({[1,1]},q);`):
fi:


if nops([args])=1 and op(1,[args])=Xavierx then
print(`Xavierx(Rset,q,x): The Generating Function for column-convex`):
print(`polyominoes (equivalenty, SAPs with at most two horizontal`):
print(`edges per vertical cross-section) according to perimeter,`):
print(` represented by q, and width of rightmost  wall, represented by x`):
print(`where interfaces are constrained to belong to the set`):
print(`Rset, for example, try: Xavier({[1,1]},q);`):
fi:

if nops([args])=1 and op(1,[args])=YafeUm then
print(`YafeUm(Umb1,F): Given an umbral operator, Umb1, in`):
print(`Maple notation, converts it to human notation`):
print(`using the symbol F for the function and D for differntiation`):
fi:

if nops([args])=1 and op(1,[args])=YafeUmSc then
print(`YafeUmSc(UmbSc,F,D): Given an umbral `):
print(`Scheme UmbSc, in`):
print(`Maple notation, converts it to human notation`):
print(`using the symbols F[1], F[2] for the functions `):
print(`and D for differntiation of respective variables `):
fi:

end:

#Alphabet(n,L,R): The set of legal paranthesis with <=n L's and n R's
Alphabet:=proc(n,L,R)
local i,gu:
gu:={}:

for i from 1 to n do
 gu:=gu union Alphabet1(i,L,R):
od:
gu:
end:



#Alphabet1(n,L,R): The set of legal paranthesis with n L's and n R's
Alphabet1:=proc(n,L,R)
local gu,k,gu1,gu2,i,j:
option remember:

if n=0 then
RETURN({[]}):
fi:

gu:={}:
for k from 0 to n-1 do
gu1:=Alphabet1(k,L,R):
gu2:=Alphabet1(n-k-1,L,R):

for i from 1 to nops(gu1) do
for j from 1 to nops(gu2) do
gu:=gu union {[L,op(gu1[i]),R,op(gu2[j])]}:
od:
od:
od:
gu:
end:


#JoinLL(LET1,ClosedI,i,L,R,C): Given a letter LET1,
# (phrased in terms of L and R and C(for closed) ), 
# A set of intervals of the form [a,a+1] indicating
# "closed blocks", and a location i, for which
# such that LET1[1][i]=LET1[1][i+1]=L does the
# operation LLpRqR->CCpLqR and appends to ClosedI 
# the interval [i,i+1]
JoinLL:=proc(LET1,ClosedI,i,L,R,C)
local w,w1,j,j1:
w:=LET1:

if not (w[i]=L and w[i+1]=L) then
ERROR(` The locations `, i , ` and `, i+1, `should be `, L):
fi:

j:=FindR(w,i,L,R,C):
j1:=FindR(w,i+1,L,R,C):

w1:=[op(1..i-1,w),C,C,op(i+2..j1-1,w),L,op(j1+1..j-1,w),op(j..nops(w),w)]:
w1,ClosedI union {[i,i+1]}:
end:


#JoinRR(LET1,ClosedI,j,L,R,C): Given a letter LET1,
# (phrased in terms of L and R and C(for closed) ), 
# A set of intervals of the form [a,a+1] indicating
# "closed blocks", and a location i, for which
# such that LET1[1][i]=LET1[1][i+1]=L does the
# operation LpLqRR->LpRqCC 
# operation LLpRqR->CCpLqR and appends to ClosedI 
# the interval [i,i+1]
JoinRR:=proc(LET1,ClosedI,j,L,R,C)
local w,w1,i,i1:
w:=LET1:

if not (w[j]=R and w[j+1]=R) then
ERROR(` The locations `, i , ` and `, i+1, `should be `, R):
fi:

i1:=FindL(w,j,L,R,C):
i:=FindL(w,j+1,L,R,C):

w1:=[op(1..i1-1,w),R,op(i1+1..j-1,w),C,C,op(j+2..nops(w),w)]:
w1,ClosedI union {[j,j+1]}:
end:


#JoinRL(LET1,ClosedI,j,L,R,C): Given a letter LET1,
# (phrased in terms of L and R), 
# such that LET1[1][j-1]=R, LET1[1][j]=L does the
# operation pRLq->pCCq and appends [j,j+1] to
# on ClosedI
JoinRL:=proc(LET1,ClosedI,j,L,R,C)
local w,w1:
w:=LET1:

if not (w[j-1]=R and w[j]=L) then
ERROR(` The locations `, j-1 , ` and `, j, `should be `, R, L):
fi:
w1:=[op(1..j-2,w),C,C,op(j+1..nops(w),w)]:

w1,ClosedI union {[j-1,j]}:
end:



#OneStepLL(LET1,ClosedI,L,R,C): Set of all pairs [LET2,ClosedI]
# that can be obtained
#from the letter LET1, with its accompanying 
#ClosedI
#L an R denote ( and ) resp. and C is "Closed"
#by performing one JoinLL
OneStepLL:=proc(LET1,ClosedI,L,R,C)
local i,w,gu:
w:=LET1:
gu:={}:
for i from 1 to nops(w)-1 do
if  w[i]=L and w[i+1]=L then
 gu:=gu union {[JoinLL(w,ClosedI,i,L,R,C)]}:
fi:
od:
gu:
end:



#OneStepRR(LET1,ClosedI,L,R,C): Set of all pairs [LET2,ClosedI]
# that can be obtained
#from the letter LET1, with its accompanying 
#ClosedI by performing one JoinRR
#L an R denote ( and ) resp. and C is "Closed"
#
OneStepRR:=proc(LET1,ClosedI,L,R,C)
local i,w,gu:
w:=LET1:

gu:={}:
for i from 1 to nops(w)-1 do
if  w[i]=R and w[i+1]=R then
 gu:=gu union {[JoinRR(w,ClosedI,i,L,R,C)]}:
fi:
od:
gu:
end:





#OneStepRL(LET1,ClosedI,L,R,C): Set of all pairs [LET2,ClosedI]
# that can be obtained
#from the letter LET1, with its accompanying 
#ClosedI by performing one JoinRL
#L an R denote ( and ) resp. and C is "Closed"
#
OneStepRL:=proc(LET1,ClosedI,L,R,C)
local i,w,gu:
w:=LET1:
gu:={}:
for i from 2 to nops(w) do
if  w[i-1]=R and w[i]=L then
 gu:=gu union {[JoinRL(w,ClosedI,i,L,R,C)]}:
fi:
od:
gu:
end:

#OneStepAll(LET1,ClosedI,L,R,C): Set of all pairs 
# [LET1, ClosedI]
# that can be obtained
#from the letter LET1, with its accompanying ClosedI
#L an R denote ( and ) resp. and C is neutral
#
OneStepAll:=proc(LET1,ClosedI,L,R,C):
OneStepLL(LET1,ClosedI,L,R,C) union
OneStepRR(LET1,ClosedI,L,R,C) union
OneStepRL(LET1,ClosedI,L,R,C):
end:




# FindR(w,i,L,R,C): Given a legal word w (in {L,R,C}), and
# a place i that is an L, find the location of its R-mate
FindR:=proc(w,i,L,R,C)
local w1,j,gu:

if w[i]<>L then
ERROR(` The `, i, `place of`, w, `should have been an `, L):
fi:

gu:=1:
w1:=subs({L=1,R=-1,C=0},w):

for j from i+1 while gu>0 do
gu:=gu+w1[j]:
od:
j-1:
end:



# FindL(w,j,L,R,C): Given a legal word w (in {L,R}), and
# a place j  that is an R, find the location of its L-mate
FindL:=proc(w,j,L,R,C)
local w1,i,gu:

if w[j]<>R then
ERROR(` The `, i, `place of`, w, `should have been an `, L):
fi:

gu:=1:
w1:=subs({L=-1,R=1,C=0},w):

for i from j-1 by -1 while gu>0 do
gu:=gu+w1[i]:
od:
i+1:
end:


#OneStepAlls(SetLet,L,R,C): Given a set of pairs
# of the form [LET1, ClosedI] finds the set of all those
# that can be obtained
#from the letters LET1, with their accompanying ClosedI
#L an R denote ( and ) resp. and C is neutral
#
OneStepAlls:=proc(SetLet,L,R,C)
local gu,i:
gu:={}:

for i from 1 to nops(SetLet) do
gu:=gu union OneStepAll(op(SetLet[i]),L,R,C):
od:
gu:
end:

#PrePreFollowers(LET,L,R,C): Given an Umbral  letter LET
#and symbols L,R,C, finds the set of all pre-pre-followers
#which are pairs of the form [LET1,ClosedI]
PrePreFollowers:=proc(LET,L,R,C)
local gu,mu:

gu:={[LET,{}]}:
mu:={[LET,{}]}:

while mu<>{} do
mu:=OneStepAlls(mu,L,R,C):
gu:=gu union mu:
od:
gu:
end:

#Hamsekh(PreLet,L,R): Given a partial pre-letter consisting
#of L,R,C where the L and R finds the three possible continuations
#of the first unconverted L or R
#with [L,-1], [L,0], [L,1] representing a Left Arm going down,
#straight, and up respectively, and similary for
#[R,-1],[R,0],[R,1]
Hamsekh:=proc(PreLet,L,R)
local mu,i:

for i from 1 to nops(PreLet) do
mu:=PreLet[i]:
if mu=L or mu=R then
RETURN({[op(1..i-1,PreLet),[mu,-1],op(i+1..nops(PreLet),PreLet)],
[op(1..i-1,PreLet),[mu,0],op(i+1..nops(PreLet),PreLet)],
[op(1..i-1,PreLet),[mu,1],op(i+1..nops(PreLet),PreLet)]}):
fi:
od:
{}:
end:

#Hamsekhs(SetPreLet,L,R): Given a set of Pre-Letters
#finds its continuations (i.e. it is the union
#of Hamsekh for every member
Hamsekhs:=proc(SetPreLet,L,R)
local gu,i:
gu:={}:

for i from 1 to nops(SetPreLet) do
gu:=gu union Hamsekh(SetPreLet[i],L,R):
od:

gu:
end:

#Hemsekhim1(LET,L,R): All the continuations of the letter LET
Hemsekhim1:=proc(LET,L,R)
local gu,mu:
mu:={LET}:

while mu<>{} do
gu:=mu:
mu:=Hamsekhs(mu,L,R):
od:
gu:
end:

#Hemsekhim(LET,L,R): All the continuations of the letter [LET,CovredI]
Hemsekhim:=proc(LET,L,R)
local mu,i:
mu:=Hemsekhim1(LET[1],L,R):
{seq([mu[i],LET[2]],i=1..nops(mu))}:
end:

#Hemsekhims(SetLET,L,R): All the continuations of the 
#members of SetLet, which are letters of the form [LET,CovredI]
Hemsekhims:=proc(SetLET,L,R)
local gu,i:
gu:={}:

for i from 1 to nops(SetLET) do
gu:=gu union Hemsekhim(SetLET[i],L,R):
od:
gu:
end:


#PreFollowers(LET,L,R,C): Given an Umbral letter in
#the SAP alphabet (a word in the alphabet {L,R})
#finds all the pre-followers
PreFollowers:=proc(LET,L,R,C):
Hemsekhims(PrePreFollowers(LET,L,R,C),L,R):
end:

#FreeSpace(PreLet):Given a pre-letter, PreLet, finds the
#list of free intervals,in order from 
#bottom to top where new pairs of LRs may be inserted
FreeSpace:=proc(PreLet)
local gu,mu,n,i,lu:
gu:=PreLet[1]:
mu:=PreLet[2]:
n:=nops(gu):
lu:=[]:

for i from 1 to n+1 do

if not member([i-1,i], mu)then 
 lu:=[op(lu),[i-1,i]]:
fi:

od:
lu:
end:

  
#Ckr(k,r):The set of vectors of non-negative
#integers of length k that sum up to r
Ckr:=proc(k,r)
local gu,mu,i,i1:
if k=0 then
 if r=0 then
  RETURN({[]}):
 else
  RETURN({}):
 fi:
fi:

gu:={}:

for i from 0 to r do
mu:=Ckr(k-1,r-i):

for i1 from 1 to nops(mu) do
 gu:=gu union {[op(mu[i1]),i]}:
od:
od:
gu:
end:



#Followers1(PreLet,C):Given a Pre-Letter, PreLet, and
#symbols L,R,C, finds all the followers 
#with up to N LR pairs
Followers1:=proc(PreLet,C,N)
local gu1,gu2,i,co,kama,FS,lu,GU,kama1:
gu1:=PreLet[1]:
gu2:=PreLet[2]:
FS:=FreeSpace(PreLet):

co:=0:

for i from 1 to nops(gu1) do
if gu1[i]=C then
 co:=co+1:
fi:
od:

kama:=N-(nops(gu1)-co)/2:

GU:={}:

for kama1 from 0 to kama do

lu:=Ckr(nops(FS),kama1):

for i from 1 to nops(lu) do
GU:=GU union {[gu1,gu2,FS,lu[i]]}:
od:

od:

GU:
end:



#Followers(LET,L,R,C,n): All the ways to follow
#the letter LET if the maximal number of pairs  of
#parantheses is n
Followers:=proc(LET,L,R,C,n)
local mu, gu,i:

mu:=PreFollowers(LET,L,R,C):
gu:={}:

for i from 1 to nops(mu) do
gu:=gu union Followers1(mu[i],C,n):
od:
gu:
end:


#OutLetter(PreLet): Given a pre-letter
#[Word,ClosedI,FS,NewLR_assignment], finds
#the resulting letter that emerges on the right
OutLetter:=proc(PreLet,L,R,C)
local gu1,i,gu3,gu4,T,n,j,w,i1:
gu1:=PreLet[1]:
gu3:=PreLet[3]:
gu4:=PreLet[4]:
n:=nops(gu1):
for i from 0 to n do
T[i]:=[]:
od:

for j from 1 to nops(gu3) do
if gu4[j]<>0 then
T[gu3[j][1]]:=[seq(op([L,R]),i1=1..gu4[j])]:
fi:
od:

w:=T[0]:

for i from 1 to n do

if gu1[i]<>C then
 w:=[op(w),gu1[i][1]]:
fi:

w:=[op(w),op(T[i])]:

od:

w:
end:


#APU(InpLet,OutPreLet,L,R,C,x,q,A): Given an input letter from the left
#InpLet, and an output Pre-letter from the right,
#OutPreLet, finds the Atomic-Pre-Umbra
#with respect to the "length-variables" x[1],x[2],...
#and q is the variable that keeps track of the length
#of the partial SAP
APU:=proc(InpLet,OutPreLet,L,R,C,x,q,A)
local w,i,gu1,gu2,gu3,gu4,mu,gu,n,T,i1,lu,T0,Tinf,xinf,lub,lut:
n:=nops(InpLet):
gu1:=OutPreLet[1]:
gu2:=OutPreLet[2]:
gu3:=OutPreLet[3]:
gu4:=OutPreLet[4]:
w:=[]:
for i from 1 to nops(gu1) do
if gu1[i]=C then
 w:=[op(w),C]:
else
 w:=[op(w),gu1[i][1]]:
fi:
od:

for i from 1 to nops(gu2) do
 T[gu2[i][1]]:=0:
od:

for i from 1 to nops(gu3) do
 T[gu3[i][1]]:=gu4[i]:
od:

gu:=[[x[0],seq(op([q*x[2*i1-1],x[2*i1]]),i1=1..T[0])]]:
mu:=2*T[0]:
for i from 1 to n do
if w[i]=C then
 gu:=[op(gu),[x[mu],seq(op([q*x[mu+2*i1+1],x[mu+2*i1+2]]),i1=0..T[i]-1)]]:
 mu:=mu+2*T[i]:
else
 gu:=[op(gu),[x[mu+1],seq(op([q*x[mu+2*i1+2],x[mu+2*i1+3]]),i1=0..T[i]-1)]]:
 mu:=mu+2*T[i]+1:
fi:
od:
gu:


for i from 1 to nops(gu) do
T[i-1]:=gu[i]:
od:

for i from 1 to nops(gu1) do
if gu1[i][1]<>C then
 if gu1[i][2]=-1 then
 T[i-1]:=[op(T[i-1]),q*T[i][1]]:
 elif gu1[i][2]=1 then
 T[i]:=[q*T[i-1][nops(T[i-1])],op(T[i])]:
 fi:
fi:
od:

gu:=[seq(T[i],i=0..nops(gu1))]:

lu:=q^(nops(OutLetter(OutPreLet,L,R,C))):

for i from 1 to nops(gu2) do
lu:=lu*q^(A[gu2[i][1]]):
od:

T0:=[op(2..nops(T[0]),T[0])]:
Tinf:=[op(1..nops(T[nops(gu1)])-1,T[nops(gu1)])]:
xinf:=op(nops(T[nops(gu1)]),T[nops(gu1)]):

for i1 from 1 to nops(T0) do
lu:=lu*T0[i1]/(1-T0[i1]):
od:

for i from 1 to nops(gu1)-1 do
 lu:=lu*PzA(T[i],A[i]):
od:

for i1 from 1 to nops(Tinf) do
lu:=lu*Tinf[i1]/(1-Tinf[i1]):
od:

lu:=normal(lu):
lu:=normal(subs({x[0]=1},lu)):
if expand(subs(xinf=1,denom(lu)))=0 then
  if expand(subs(xinf=1,numer(lu)))<>0 then
   ERROR(`Something is wrong`):
  fi:

 lut:= subs(xinf=1,diff(numer(lu),xinf)):
 lub:= subs(xinf=1,diff(denom(lu),xinf)):
 lu:=normal(lut/lub):
else
lu:=normal(subs({xinf=1},lu)):
fi:
lu:

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)),
symbolic)):
end:

 

#PreUmbra(LETTER,L,R,n,x,q,A,Z,HATKH,SOF)
#Given a letter LETTER, in terms of legal
#parantheses using L and R
#and an integer n specifying the max. number of
#LR-pairs, and a letter Z
#and a symbol A, finds the
#pre-umbra (i.e. the image of a generic monomial
#x[1]^A[1]*x[2]^A[2]*...*x[r]^A[r] where r is nops(LETTER)-1
#as a linear combination of Z[LETTER1], where LETTER1
#ranges over all followers of LETTER

PreUmbra:=proc(LETTER,L,R,n,x,q,A,Z,HATKH,SOF)
local gu,mu,C,i,n1,lu,tu,i1:

if LETTER=SOF then
RETURN(0):
fi:

if LETTER=HATKH then
 gu:=0:
 for i from 1 to n do
  lu:=[seq(op([L,R]),i1=1..i)]:
  tu:=q^3*x[1]/(1-q*x[1]):
   for i1 from 1 to i-1 do
    tu:=tu*x[2*i1]/(1-x[2*i1])*q^3*x[2*i1+1]/(1-q*x[2*i1+1]):
   od:
   gu:=gu+tu*Z[lu]:
 od:
RETURN(gu):
fi:
  
mu:=Followers(LETTER,L,R,C,n):

gu:=0:

for i from 1 to nops(mu) do
gu:=gu+APU(LETTER,mu[i],L,R,C,x,q,A)*Z[OutLetter(mu[i],L,R,C)]:
od:

n1:=nops(LETTER)/2:
if LETTER=[L,seq(op([L,R]),i=1..n1-1),R] then
lu:=0:
 for i from 1 to n1 do
  lu:=lu+A[2*i-1]:
 od:

  gu:=gu+Z[SOF]*q^lu:
fi:

gu:
end:

#SAPUmSc(n,x,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
#followed by the initial values of the corresponding function
#followed by the set of final letters (or states)
#for self-avoiding polygons with the property
#that each vertical cross section has at most n LR-pairs
SAPUmSc:=proc(n,x,q)
local mu,L,R,A,Z,HATKH,SOF,i1,S,T,U,alist,xlist,gu,gu1,i,j,lu,tu,m,mu1,
yu,zu1:
mu:=Alphabet(n,L,R) union {HATKH,SOF}:
gu:=[]:

for m from 1 to n do

mu1:=Alphabet1(m,L,R):
for i from 1 to nops(mu1) do
xlist:=[seq(x[i1],i1=1..2*m-1)]:
alist:=[seq(A[i1],i1=1..2*m-1)]:
S[mu1[i]]:=xlist:
 U[mu1[i]]:=0:

lu:=PreUmbra(mu1[i],L,R,n,x,q,A,Z,HATKH,SOF):

for j from 1 to nops(mu) do
T[mu1[i],mu[j]]:=ToUmbra(coeff(lu,Z[mu[j]],1),xlist,alist):
od:


od:

od:

mu1:={HATKH,SOF}:
for i from 1 to nops(mu1) do
xlist:=[]:
alist:=[]:
S[mu1[i]]:=xlist:

if mu1[i]=HATKH then
 U[mu1[i]]:=1:
else
 U[mu1[i]]:=0:
fi:

lu:=PreUmbra(mu1[i],L,R,n,x,q,A,Z,HATKH,SOF):

for j from 1 to nops(mu) do
T[mu1[i],mu[j]]:=ToUmbra(coeff(lu,Z[mu[j]],1),xlist,alist):
od:
od:


gu:=[]:
tu:=[]:

mu:=convert(mu,list):

yu:=[]:

for i from 1 to nops(mu) do
 tu:=[op(tu),S[mu[i]]]:
 yu:=[op(yu),U[mu[i]]]:
 if mu[i]=SOF then
  zu1:={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:

#############Programs Taken 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(SimplifyUmbra({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), symbolic):
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,pu:
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):
pu:=normal(coeff(gu,F(evar))):
if pu<>0 then
mu:=mu union {[factor(pu),evar[1],evar[2]]}:
fi:

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:






# 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:




#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:

#############End Programs Taken from ROTA###############


# 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(normal(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,normal(mu)):

Sid:=[op(Sid),normal(mu)]:

if i mod 2=0 then
lprint([seq(Sid[2+2*i1],i1=1..i/2-1)]):
fi:

od:

lu,[seq(Sid[2+2*i1],i1=1..n/2-1)]:

end:

#Milim1(L,R,C,n,M): All possible templates of self-avoiding
#polygons with at most n LR's and of length M
Milim1:=proc(L,R,C,n,M)
local gu,lu,i,j,mu1,gu1,mu:
option remember:

if M=1 then
lu:=[L,R]:
gu:={[lu]}:

for i from 2 to n do
 lu:=[op(lu),L,R]:
 gu:=gu union {[lu]}:
od:
RETURN(gu):
fi:

mu:=Milim1(L,R,C,n,M-1):

gu:={}:

for i from 1 to nops(mu) do
 mu1:=mu[i]:
 gu1:=Followers(mu1[nops(mu1)],L,R,C,n):

for j from 1 to nops(gu1) do
 gu:=gu union {[op(mu1),gu1[j],OutLetter(gu1[j],L,R,C)]}:
od:
od:
gu:

end:



#Milim(L,R,C,n,M): all terminal words in the SAP alphabet
#of length M and <=n LR's 
Milim:=proc(L,R,C,n,M)
local mu,gu,i,n1,i1,mu1,akha:
mu:=Milim1(L,R,C,n,M):
gu:={}:

for i from 1 to nops(mu) do
mu1:=mu[i]:
akha:=mu1[nops(mu1)]:
n1:=nops(akha)/2:
 if akha=[L,seq(op([L,R]),i1=1..n1-1),R] then
   gu:=gu union {mu1}:
 fi:
od:
gu:
end:


SAPUS1:=
[{3}, [[{}, {}, {}], [{[-q**3*x[1]/(-1+q*x[1]), [], []]}, {[q**2/(-x[1]+q)*x[
1], [1], [q]], [q**2*x[1]**2*(q-1)**2*(q+1)**2/(-x[1]+q)**2/(-1+q*x[1])**2, [0]
, [x[1]]], [-q**2*x[1]**2*(-q*x[1]+2*q**2-1)/(-x[1]+q)**2/(-1+q*x[1]), [0], [q]
]}, {}], [{}, {[1, [0], [q]]}, {}]], [1, 0, 0], [[], [x[1]], []]]:

SAPUS2:=
[{5}, [[{}, {}, {}, {}, {}], [{[-q**3*x[1]/(-1+q*x[1]), [], []]}, {[q**2/(-x[1]
+q)*x[1], [1], [q]], [q**2*x[1]**2*(q-1)**2*(q+1)**2/(-x[1]+q)**2/(-1+q*x[1])**
2, [0], [x[1]]], [-q**2*x[1]**2*(-q*x[1]+2*q**2-1)/(-x[1]+q)**2/(-1+q*x[1]), [0
], [q]]}, {[-q**2*x[1]**2*(q-1)*(q+1)/(-x[1]+q)**2/(-1+q*x[1]), [0, 0, 0], [x[1
], q*x[1], q]], [q**2*x[1]**2*(q-1)**2*(q+1)**2/(-x[1]+q)**2/(-1+q*x[1])**2, [0
, 0, 0], [x[1], q*x[1], x[1]]], [-q**2*x[1]**2*(q-1)*(q+1)/(-x[1]+q)**2/(-1+q*x
[1]), [0, 0, 0], [q, q*x[1], x[1]]], [q**2/(-x[1]+q)**2*x[1]**2, [0, 0, 0], [q
, q*x[1], q]]}, {[-2*q**2*x[1]**2*(-q*x[1]+2*q**2-1)/(-x[1]+q)**2/(-1+q*x[1]), 
[0, 0, 0], [q, 1, q]], [q**2*x[1]**2*(q-1)**2*(q+1)**2/(-x[1]+q)**2/(-1+q*x[1])
**2, [0, 0, 0], [x[1], 1, q]], [-q**2*x[1]*(q-1)*(q+1)/(-x[1]+q)/(-1+q*x[1])**2
, [0, 0, 0], [x[1], q*x[1], q]], [q**2*x[1]**2*(q-1)**2*(q+1)**2/(-x[1]+q)**2/(
-1+q*x[1])**2, [0, 0, 0], [q, 1, x[1]]], [-q**2*x[1]*(q-1)*(q+1)/(-x[1]+q)/(-1+
q*x[1])**2, [0, 0, 0], [q, q*x[1], x[1]]], [q**2/(-x[1]+q)*x[1], [0, 0, 1], [q
, 1, q]], [q**2/(-x[1]+q)*x[1], [1, 0, 0], [q, 1, q]], [2*q**2/(-x[1]+q)/(-1+q*
x[1])*x[1], [0, 0, 0], [q, q*x[1], q]]}, {}], [{[-q**6*x[1]*x[2]*x[3]/(-1+q*x[1
])/(-1+x[2])/(-1+q*x[3]), [], []]}, {[(q-1)*(q+1)*(q**2*x[2]-1)*x[1]**2*q**5*x[
2]*x[3]/(-1+q*x[1])**2/(-x[1]+q*x[2])/(-x[1]+q)/(-1+x[2])/(-1+q*x[3]), [0], [x[
1]]], [-x[1]*x[3]*x[2]*q**5*(x[3]*q**2*x[2]-2*q*x[1]*x[3]+x[3]*x[2]-2*q*x[2]**2
+q**2*x[2]*x[1]+x[2]*x[1])/(-q*x[2]+x[3])/(-1+q*x[1])/(-1+q*x[3])/(-1+x[2])**2/
(-x[1]+q*x[2]), [0], [q*x[2]]], [x[1]*x[3]*x[2]*q**5*(x[3]*q**2*x[2]-2*x[1]*x[3
]*q*x[2]+x[3]+q**2*x[2]*x[1]+x[1]-2*q)/(x[3]-q)/(-1+q*x[1])/(-1+q*x[3])/(-1+x[2
])**2/(-x[1]+q), [0], [q]], [(q-1)*(q+1)*(q**2*x[2]-1)*q**5*x[1]*x[2]*x[3]**2/(
-1+q*x[3])**2/(-q*x[2]+x[3])/(-1+q*x[1])/(-1+x[2])/(x[3]-q), [0], [x[3]]]}, {[-
q**5*x[1]**2*x[2]*x[3]*(q**2*x[2]-1)/(-x[1]+q)/(-x[1]+q*x[2])/(-1+x[2])/(-1+q*x
[3])/(-1+q*x[1]), [0, 0, 0], [q, q*x[1], x[1]]], [-q**5*x[1]*x[2]*x[3]**2*(q**2
*x[2]-1)/(-1+q*x[1])/(x[3]-q)/(-q*x[2]+x[3])/(-1+x[2])/(-1+q*x[3]), [0, 0, 0], 
[x[3], q*x[3], q]], [(q-1)*(q+1)*(q**2*x[2]-1)*x[1]**2*q**5*x[2]*x[3]/(-1+q*x[1
])**2/(-x[1]+q*x[2])/(-x[1]+q)/(-1+x[2])/(-1+q*x[3]), [0, 0, 0], [x[1], q*x[1]
, x[1]]], [-(-1+q*x[1])*(-1+q*x[3])*q**4*x[2]**2*x[1]*x[3]/(q*x[3]-x[2])/(-1+x[
2])**2/(-x[1]+q)/(q*x[1]-x[2])/(x[3]-q), [0, 0, 0], [q, x[2], q]], [x[3]*x[2]*x
[1]*q**4*(q*x[1]-x[2])*(-1+q*x[3])/(-1+q*x[1])/(x[3]-q)/(q*x[3]-x[2])/(-1+x[2])
**2/(-x[1]+q*x[2]), [0, 0, 0], [q*x[2], x[2], q]], [-x[3]*x[2]*x[1]*q**4*(q-1)*
(q+1)*(-x[2]*x[1]+x[3])/(-1+x[2])/(-1+q*x[3])/(x[3]-q)/(-x[1]+x[3])/(-x[1]+q)/(
q*x[3]-x[2]), [0, 0, 0], [q, q*x[3], x[3]]], [x[3]*x[2]*x[1]*q**4*(-1+q*x[1])*(
q*x[3]-x[2])/(q*x[1]-x[2])/(-1+q*x[3])/(-q*x[2]+x[3])/(-1+x[2])**2/(-x[1]+q), [
0, 0, 0], [q, x[2], q*x[2]]], [x[2]*x[1]*x[3]**2*q**4*(q-1)*(q+1)*(-x[1]+q*x[2]
)/(-q*x[2]+x[3])/(x[3]-q)/(-1+q*x[3])/(-x[1]+x[3])/(q*x[1]-x[2])/(-x[1]+q), [0
, 0, 0], [q, q*x[1], x[3]]], [(x[3]**2*q**3*x[2]-q**2*x[2]**2*x[3]+q**2*x[2]*x[
1]-q**2*x[3]**2*x[1]-q*x[1]*x[3]-q*x[2]**2+x[1]*x[3]*q*x[2]+x[3]*x[2])*q**4*x[1
]*x[2]*x[3]/(-q*x[2]+x[3])/(-x[1]+q*x[2])/(-1+q*x[1])/(x[3]-q)/(-1+x[2])/(q*x[3
]-x[2]), [0, 0, 0], [q*x[2], q*x[3], q]], [x[3]*x[2]*x[1]*q**4*(x[3]*x[2]-x[1])
/(-x[1]+x[3])/(-1+x[2])/(-x[1]+q)/(q*x[1]-x[2])/(x[3]-q), [0, 0, 0], [q, q*x[1]
, q]], [(-q**3*x[1]**2*x[2]-x[3]*q**2*x[2]+q**2*x[1]**2*x[3]+q**2*x[2]**2*x[1]+
q*x[1]*x[3]+q*x[2]**2-x[1]*x[3]*q*x[2]-x[2]*x[1])*q**4*x[1]*x[2]*x[3]/(-1+q*x[3
])/(-x[1]+q*x[2])/(-1+x[2])/(-x[1]+q)/(q*x[1]-x[2])/(-q*x[2]+x[3]), [0, 0, 0], 
[q, q*x[1], q*x[2]]], [x[3]*x[2]*x[1]*q**4*(-x[2]*x[1]+x[3])/(-1+x[2])/(x[3]-q)
/(-x[1]+x[3])/(q*x[3]-x[2])/(-x[1]+q), [0, 0, 0], [q, q*x[3], q]], [-x[2]*x[1]*
x[3]**2*q**4*(q-1)**2*(q+1)**2*(-x[1]+q*x[2])/(-x[1]+x[3])/(-q*x[2]+x[3])/(x[3]
-q)/(-x[1]+q)/(-1+q*x[1])/(q*x[1]-x[2])/(-1+q*x[3]), [0, 0, 0], [x[1], q*x[1], 
x[3]]], [-(q-1)*(q+1)*(x[3]**2*q**3*x[2]-q**2*x[2]**2*x[3]+q**2*x[2]*x[1]-q**2*
x[3]**2*x[1]-q*x[1]*x[3]-q*x[2]**2+x[1]*x[3]*q*x[2]+x[3]*x[2])*q**4*x[1]*x[2]*x
[3]/(-q*x[2]+x[3])/(-1+q*x[1])/(-1+q*x[3])/(x[3]-q)/(q*x[3]-x[2])/(-1+x[2])/(-x
[1]+q*x[2]), [0, 0, 0], [q*x[2], q*x[3], x[3]]], [-x[3]*x[1]**2*x[2]*q**4*(q-1)
*(q+1)*(-q*x[2]+x[3])/(-1+q*x[1])/(-x[1]+x[3])/(x[3]-q)/(q*x[3]-x[2])/(-x[1]+q*
x[2])/(-x[1]+q), [0, 0, 0], [x[1], q*x[3], q]], [-(q-1)**2*(q+1)**2*(-1+q*x[1])
*q**4*x[1]*x[3]**2*x[2]**2/(x[3]-q)/(q*x[3]-x[2])/(-1+x[2])/(-x[1]+q)/(q*x[1]-x
[2])/(-q*x[2]+x[3])/(-1+q*x[3]), [0, 0, 0], [q, x[2], x[3]]], [(q-1)*(q+1)*(q**
2*x[2]-1)*q**5*x[1]*x[2]*x[3]**2/(-1+q*x[3])**2/(-q*x[2]+x[3])/(-1+q*x[1])/(-1+
x[2])/(x[3]-q), [0, 0, 0], [x[3], q*x[3], x[3]]], [-x[3]*x[2]*x[1]*q**4*(q-1)*(
q+1)*(x[3]*x[2]-x[1])/(-1+x[2])/(-x[1]+q)/(x[3]-q)/(-x[1]+x[3])/(q*x[1]-x[2])/(
-1+q*x[1]), [0, 0, 0], [x[1], q*x[1], q]], [-(q*x[1]-x[2])*(q*x[3]-x[2])*q**4*x
[1]*x[3]/(-1+q*x[3])/(-1+q*x[1])/(-x[1]+q*x[2])/(-1+x[2])**2/(-q*x[2]+x[3]), [0
, 0, 0], [q*x[2], x[2], q*x[2]]], [-x[3]*x[1]**2*x[2]*q**4*(q-1)**2*(q+1)**2*(q
*x[3]-x[2])/(-x[1]+q)/(-1+q*x[3])/(-x[1]+q*x[2])/(-1+x[2])/(-1+q*x[1])/(q*x[1]-
x[2])/(-q*x[2]+x[3]), [0, 0, 0], [x[1], x[2], q*x[2]]], [x[2]*x[1]*x[3]**2*q**4
*(q-1)**2*(q+1)**2*(q*x[1]-x[2])/(q*x[3]-x[2])/(-1+x[2])/(x[3]-q)/(-1+q*x[1])/(
-x[1]+q*x[2])/(-q*x[2]+x[3])/(-1+q*x[3]), [0, 0, 0], [q*x[2], x[2], x[3]]], [x[
2]**2*x[1]**2*x[3]**2*q**4*(q-1)**4*(q+1)**4/(-1+q*x[3])/(-1+q*x[1])/(-x[1]+q*x
[2])/(-q*x[2]+x[3])/(x[3]-q)/(q*x[1]-x[2])/(-x[1]+q)/(q*x[3]-x[2]), [0, 0, 0], 
[x[1], x[2], x[3]]], [-(q-1)*(q+1)*(-q**3*x[1]**2*x[2]-x[3]*q**2*x[2]+q**2*x[1]
**2*x[3]+q**2*x[2]**2*x[1]+q*x[1]*x[3]+q*x[2]**2-x[1]*x[3]*q*x[2]-x[2]*x[1])*q
**4*x[1]*x[2]*x[3]/(-x[1]+q)/(-1+q*x[3])/(-x[1]+q*x[2])/(-1+x[2])/(-1+q*x[1])/(
q*x[1]-x[2])/(-q*x[2]+x[3]), [0, 0, 0], [x[1], q*x[1], q*x[2]]], [x[3]*x[1]**2*
x[2]*q**4*(q-1)**2*(q+1)**2*(-q*x[2]+x[3])/(x[3]-q)/(-1+q*x[3])/(-x[1]+x[3])/(-
1+q*x[1])/(-x[1]+q*x[2])/(-x[1]+q)/(q*x[3]-x[2]), [0, 0, 0], [x[1], q*x[3], x[3
]]], [(q-1)**2*(q+1)**2*(-1+q*x[3])*q**4*x[2]**2*x[1]**2*x[3]/(-1+q*x[1])/(q*x[
3]-x[2])/(q*x[1]-x[2])/(x[3]-q)/(-1+x[2])/(-x[1]+q)/(-x[1]+q*x[2]), [0, 0, 0], 
[x[1], x[2], q]]}, {[-x[1]*x[3]*x[2]*q**5*(x[3]*q**2*x[2]-2*q*x[1]*x[3]+x[3]*x[
2]-2*q*x[2]**2+q**2*x[2]*x[1]+x[2]*x[1])/(-q*x[2]+x[3])/(-1+q*x[1])/(-1+q*x[3])
/(-1+x[2])**2/(-x[1]+q*x[2]), [0, 0, 0], [q*x[2], x[2], q*x[2]]], [x[3]*x[2]*q
**4*(-x[3]*x[1]**2*q-q**3*x[2]*x[3]-x[3]*q**2*x[1]**3*x[2]+2*x[1]**2*x[2]*q*x[3
]-q**2*x[1]*x[2]**2*x[3]+x[3]*x[2]**2*x[1]**2*q**3+q**2*x[3]*x[1]+q**3*x[1]**2*
x[3]-x[3]*q**2*x[1]**3+x[2]*x[1]**3*q**3-x[2]**3*q**4*x[1]**2+x[2]**2*q**3*x[1]
**3-q*x[1]*x[2]**2-q**4*x[2]*x[1]**2-x[2]*x[1]*q+q**2*x[2]**2+x[1]**2*x[2]+q**3
*x[1]*x[2]**2+q**3*x[2]**3*x[1]-2*q**2*x[2]**2*x[1]**2)/(-x[1]+q)/(-1+q*x[3])/(
-x[1]+q*x[2])/(q*x[1]-x[2])/(-1+q*x[1])/(-1+x[2])/(-q*x[2]+x[3]), [0, 0, 0], [q
, q*x[1], q*x[2]]], [-x[1]*q**5*x[3]**2*x[2]*(q-1)**2*(q+1)**2/(-q*x[2]+x[3])/(
x[3]-q)/(-1+q*x[3])/(-1+x[2])/(q*x[1]-x[2])/(q*x[3]-x[2]), [0, 0, 0], [q, x[2]
, x[3]]], [(q-1)*(q+1)*(q**2*x[2]-1)*x[1]**2*q**5*x[2]*x[3]/(-1+q*x[1])**2/(-x[
1]+q*x[2])/(-x[1]+q)/(-1+x[2])/(-1+q*x[3]), [0, 0, 0], [x[1], 1, q]], [(q-1)*(q
+1)*(-q*x[2]+x[3])*q**4*x[2]*x[1]**2/(-x[1]+q*x[2])/(-1+q*x[1])/(-x[1]+x[3])/(q
*x[3]-x[2])/(-1+q*x[3])/(-x[1]+q), [0, 0, 0], [x[1], q*x[3], q]], [-x[3]*x[2]*x
[1]**2*q**5*(q-1)**2*(q+1)**2/(-x[1]+q*x[2])/(-1+q*x[1])/(q*x[1]-x[2])/(-1+x[2]
)/(q*x[3]-x[2])/(-x[1]+q), [0, 0, 0], [x[1], x[2], q]], [(q-1)*(q+1)*q**4*x[1]*
x[2]*x[3]/(-1+q*x[1])**2/(-x[1]+x[3])/(q*x[1]-x[2])/(-x[1]+q), [0, 0, 0], [x[1]
, q*x[1], q]], [2*x[1]*x[3]*x[2]*q**5*(x[3]*q**2*x[2]-2*x[1]*x[3]*q*x[2]+x[3]+q
**2*x[2]*x[1]+x[1]-2*q)/(x[3]-q)/(-1+q*x[1])/(-1+q*x[3])/(-1+x[2])**2/(-x[1]+q)
, [0, 0, 0], [q, 1, q]], [(q-1)**2*(q+1)**2*q**5*x[1]*x[2]**2*x[3]**2/(q*x[3]-x
[2])/(x[3]-q)/(-1+q*x[1])/(-1+x[2])/(-q*x[2]+x[3])/(-1+q*x[3]), [0, 0, 0], [q*x
[2], x[2], x[3]]], [(q-1)**2*(q+1)**2*x[1]**2*x[3]*q**5*x[2]**2/(-x[1]+q*x[2])/
(-x[1]+q)/(q*x[1]-x[2])/(-1+q*x[1])/(-1+x[2])/(-1+q*x[3]), [0, 0, 0], [x[1], x[
2], q*x[2]]], [-q**5*x[3]*x[2]*x[1]*(q**2*x[2]-1)/(-1+x[2])/(-1+q*x[3])/(-1+q*x
[1])**2/(-x[1]+q*x[2]), [0, 0, 0], [q, q*x[1], x[1]]], [(q**2*x[1]+q**2*x[2]**3
*x[3]-q*x[1]*x[3]-q*x[2]-q*x[2]**3-x[3]*x[2]**2*x[1]*q+x[2]**2*x[1]+x[3]*x[2])*
q**5*x[1]*x[3]/(x[3]-q)/(-1+x[2])**2/(q*x[3]-x[2])/(-1+q*x[1])/(-x[1]+q*x[2]), 
[0, 0, 0], [q*x[2], x[2], q]], [(q**2*x[1]*x[2]**3+q**2*x[3]-q*x[2]**3-x[3]*x[2
]**2*x[1]*q-q*x[2]-q*x[1]*x[3]+x[2]**2*x[3]+x[2]*x[1])*q**5*x[1]*x[3]/(-x[1]+q)
/(-1+q*x[3])/(-q*x[2]+x[3])/(-1+x[2])**2/(q*x[1]-x[2]), [0, 0, 0], [q, x[2], q*
x[2]]], [-(q**2*x[3]-2*q*x[1]*x[3]+x[3]-2*q+q**2*x[1]+x[1])*q**5*x[1]*x[2]*x[3]
/(x[3]-q)/(q*x[1]-x[2])/(q*x[3]-x[2])/(-1+x[2])**2/(-x[1]+q), [0, 0, 0], [q, x[
2], q]], [(q-1)*(q+1)*(q**2*x[2]-1)*x[1]**2*q**5*x[2]*x[3]/(-1+q*x[1])**2/(-x[1
]+q*x[2])/(-x[1]+q)/(-1+x[2])/(-1+q*x[3]), [0, 0, 0], [q, 1, x[1]]], [q**5*x[1]
*x[2]**2*x[3]*(q-1)*(q+1)/(-1+q*x[1])/(-1+x[2])/(q*x[3]-x[2])/(-1+q*x[3])/(x[3]
-q), [0, 0, 0], [q*x[2], q*x[3], x[3]]], [-x[1]*x[2]*q**4*(-q**4*x[3]**2*x[2]**
3+x[2]**3*q**3*x[3]-2*q**2*x[3]**2*x[2]**2+q**2*x[2]**2+x[2]**2*q**3*x[3]**3+q
**3*x[3]**2*x[2]**2*x[1]-q**2*x[1]*x[2]**2*x[3]-x[3]*x[2]**2*q+q**3*x[2]**2*x[3
]-x[2]*q**2*x[3]**3*x[1]-q**4*x[3]**2*x[2]-q**3*x[2]*x[1]+2*q*x[3]**2*x[2]*x[1]
+x[3]**2*x[2]+q**3*x[3]**3*x[2]-x[3]*q*x[2]+q**3*x[3]**2*x[1]+q**2*x[3]*x[1]-x[
3]**2*x[1]*q-q**2*x[3]**3*x[1])/(-q*x[2]+x[3])/(x[3]-q)/(-1+q*x[1])/(-x[1]+q*x[
2])/(-1+q*x[3])/(q*x[3]-x[2])/(-1+x[2]), [0, 0, 0], [q*x[2], q*x[3], q]], [x[2]
*x[1]*q**5*x[3]*(q**2*x[2]-1)/(-1+q*x[1])/(-1+x[2])/(-q*x[2]+x[3])/(-1+q*x[3])
**2, [0, 0, 0], [x[3], q*x[3], q]], [-x[1]*x[3]*x[2]*q**5*(x[3]*q**2*x[2]-2*q*x
[1]*x[3]+x[3]*x[2]-2*q*x[2]**2+q**2*x[2]*x[1]+x[2]*x[1])/(-q*x[2]+x[3])/(-1+q*x
[1])/(-1+q*x[3])/(-1+x[2])**2/(-x[1]+q*x[2]), [0, 0, 0], [q, 1, q*x[2]]], [-q**
5*x[1]*x[2]**2*x[3]*(q-1)*(q+1)/(-x[1]+q)/(q*x[1]-x[2])/(-1+q*x[1])/(-1+x[2])/(
-1+q*x[3]), [0, 0, 0], [x[1], q*x[1], q*x[2]]], [x[2]*x[3]*q**4*(x[3]*q*x[2]+x[
1]*x[3]-2*x[3]*x[2]*x[1]-2*q*x[1]+x[1]**2+x[2]*x[1]*q)/(x[3]-q)/(-1+q*x[1])/(-x
[1]+x[3])/(q*x[1]-x[2])/(-1+x[2])/(-x[1]+q), [0, 0, 0], [q, q*x[1], q]], [(q-1)
*(q+1)*(-x[1]+q*x[2])*q**4*x[3]**2*x[2]/(-q*x[2]+x[3])/(x[3]-q)/(-1+q*x[3])/(-1
+q*x[1])/(-x[1]+x[3])/(q*x[1]-x[2]), [0, 0, 0], [q, q*x[1], x[3]]], [-x[1]*x[3]
*x[2]*q**5*(x[3]*q**2*x[2]-2*q*x[1]*x[3]+x[3]*x[2]-2*q*x[2]**2+q**2*x[2]*x[1]+x
[2]*x[1])/(-q*x[2]+x[3])/(-1+q*x[1])/(-1+q*x[3])/(-1+x[2])**2/(-x[1]+q*x[2]), [
0, 0, 0], [q*x[2], 1, q]], [(q-1)*(q+1)*(q**2*x[2]-1)*q**5*x[1]*x[2]*x[3]**2/(-
1+q*x[3])**2/(-q*x[2]+x[3])/(-1+q*x[1])/(-1+x[2])/(x[3]-q), [0, 0, 0], [x[3], 1
, q]], [(q-1)*(q+1)*q**4*x[1]*x[2]*x[3]/(x[3]-q)/(-x[1]+x[3])/(q*x[3]-x[2])/(-1
+q*x[3])**2, [0, 0, 0], [q, q*x[3], x[3]]], [(q-1)*(q+1)*(q**2*x[2]-1)*q**5*x[1
]*x[2]*x[3]**2/(-1+q*x[3])**2/(-q*x[2]+x[3])/(-1+q*x[1])/(-1+x[2])/(x[3]-q), [0
, 0, 0], [q, 1, x[3]]], [-x[1]*x[2]*q**4*(x[2]*x[1]*q+x[3]*q*x[2]-2*q*x[3]-2*x[
3]*x[2]*x[1]+x[1]*x[3]+x[3]**2)/(x[3]-q)/(-x[1]+x[3])/(-x[1]+q)/(-1+q*x[3])/(q*
x[3]-x[2])/(-1+x[2]), [0, 0, 0], [q, q*x[3], q]]}, {}], [{}, {[q**4*x[1]*x[2]*x
[3]/(-x[1]+q)/(-1+x[2])/(x[3]-q), [1], [q]], [-(-q*x[3]-q*x[1]+2*q**2-q**2*x[2]
-3*q**3*x[1]-q**2*x[2]*x[1]**2+4*q**2*x[3]*x[1]-2*x[3]*x[2]*q**2*x[1]+2*q**2*x[
1]**2-3*q**4*x[1]*x[2]*x[3]+2*q**3*x[2]*x[1]-x[3]**2*x[1]*q-x[3]**2*x[2]*q**2+x
[3]*x[2]*x[1]-x[3]*x[1]**2*q+2*x[3]**2*x[2]*q**3*x[1]-3*q**3*x[3]**2*x[1]+2*x[3
]**2*x[1]**2*q**2+2*q**3*x[2]*x[3]+2*q**2*x[3]**2-3*q**3*x[3]+2*q**3*x[2]*x[1]
**2*x[3]-x[3]**2*x[2]*x[1]**2*q**2-3*q**3*x[1]**2*x[3]+4*q**4*x[1]*x[3])*q**4*x
[1]*x[2]*x[3]/(-1+q*x[3])/(-1+q*x[1])/(-x[1]+q)**2/(-1+x[2])**2/(x[3]-q)**2, [0
], [q]], [-(q*x[1]-x[2])*(q*x[3]-x[2])*q**4*x[1]*x[3]/(-1+q*x[3])/(-1+q*x[1])/(
-x[1]+q*x[2])/(-1+x[2])**2/(-q*x[2]+x[3]), [0], [q*x[2]]], [(q-1)*(q+1)*(q**2*x
[3]-x[1])*q**5*x[1]*x[3]*x[2]/(-1+q*x[3])/(-1+q*x[1])/(-x[1]+q)**2/(-x[1]+x[3])
/(-x[1]+q*x[2]), [0], [x[1]]], [-(q-1)*(q+1)*(-q**2*x[1]+x[3])*q**5*x[1]*x[3]*x
[2]/(-1+q*x[3])/(-1+q*x[1])/(-x[1]+x[3])/(-q*x[2]+x[3])/(x[3]-q)**2, [0], [x[3]
]]}, {[-x[1]*x[3]**2*q**4*x[2]/(x[3]-q)**2/(-1+x[2])/(-x[1]+q), [0, 0, 0], [q, 
q*x[3], q]], [(q*x[1]-x[2])*q**4*x[3]**2*x[1]/(-1+x[2])/(-1+q*x[1])/(x[3]-q)/(-
x[1]+q*x[2])/(-q*x[2]+x[3]), [0, 0, 0], [q*x[2], q*x[3], q]], [(q-1)*(q+1)*x[1]
*x[3]**2*q**4*x[2]/(-1+q*x[3])/(x[3]-q)**2/(-1+x[2])/(-x[1]+q), [0, 0, 0], [q, 
q*x[3], x[3]]], [-(q-1)*(q+1)*(-q**2*x[1]+x[3])*q**5*x[1]*x[3]*x[2]/(-1+q*x[3])
/(-1+q*x[1])/(-x[1]+x[3])/(-q*x[2]+x[3])/(x[3]-q)**2, [0, 0, 0], [x[3], q*x[3]
, x[3]]], [(q-1)*(q+1)*(q**2*x[3]-x[1])*q**5*x[1]*x[3]*x[2]/(-1+q*x[3])/(-1+q*x
[1])/(-x[1]+q)**2/(-x[1]+x[3])/(-x[1]+q*x[2]), [0, 0, 0], [x[1], q*x[1], x[1]]]
, [(-q**2*x[1]+x[3])*q**5*x[2]*x[1]*x[3]/(-1+q*x[1])/(x[3]-q)**2/(-x[1]+x[3])/(
-q*x[2]+x[3]), [0, 0, 0], [x[3], q*x[3], q]], [q**4*x[1]**2*x[2]*x[3]/(-1+x[2])
/(x[3]-q)/(-x[1]+q)**2, [0, 0, 0], [q, q*x[1], q]], [-(q-1)**2*(q+1)**2*q**5*x[
1]*x[3]**2*x[2]/(-x[1]+q)/(x[3]-q)/(-1+q*x[3])/(-x[1]+x[3])/(-x[1]+q*x[2])/(-1+
q*x[1]), [0, 0, 0], [x[1], q*x[3], x[3]]], [-(q-1)*(q+1)*(q*x[1]-x[2])*q**4*x[1
]*x[3]**2/(-1+q*x[3])/(-1+q*x[1])/(-x[1]+q*x[2])/(-1+x[2])/(-q*x[2]+x[3])/(x[3]
-q), [0, 0, 0], [q*x[2], q*x[3], x[3]]], [(q-1)*(q+1)*q**5*x[2]*x[1]**2*x[3]/(x
[3]-q)/(-1+q*x[3])/(-x[1]+x[3])/(-q*x[2]+x[3])/(-x[1]+q), [0, 0, 0], [q, q*x[1]
, x[3]]], [-(q-1)**2*(q+1)**2*q**5*x[2]*x[1]**2*x[3]/(-1+q*x[3])/(-x[1]+q)/(-x[
1]+x[3])/(-q*x[2]+x[3])/(x[3]-q)/(-1+q*x[1]), [0, 0, 0], [x[1], q*x[1], x[3]]]
, [(q-1)*(q+1)*(q*x[3]-x[2])*q**4*x[3]*x[1]**2/(-1+q*x[3])/(-x[1]+q)/(-x[1]+q*x
[2])/(-1+x[2])/(-q*x[2]+x[3])/(-1+q*x[1]), [0, 0, 0], [x[1], q*x[1], q*x[2]]], 
[-(q**2*x[3]-x[1])*q**5*x[2]*x[1]*x[3]/(-1+q*x[3])/(-x[1]+x[3])/(-x[1]+q*x[2])/
(-x[1]+q)**2, [0, 0, 0], [q, q*x[1], x[1]]], [(q-1)*(q+1)*q**5*x[1]*x[3]**2*x[2
]/(-x[1]+q)/(-1+q*x[1])/(x[3]-q)/(-x[1]+x[3])/(-x[1]+q*x[2]), [0, 0, 0], [x[1]
, q*x[3], q]], [-(q-1)*(q+1)*q**4*x[1]**2*x[2]*x[3]/(-1+q*x[1])/(-1+x[2])/(x[3]
-q)/(-x[1]+q)**2, [0, 0, 0], [x[1], q*x[1], q]], [-(q*x[3]-x[2])*q**4*x[3]*x[1]
**2/(-q*x[2]+x[3])/(-1+x[2])/(-x[1]+q*x[2])/(-x[1]+q)/(-1+q*x[3]), [0, 0, 0], [
q, q*x[1], q*x[2]]]}, {[(q-1)*(q+1)*(q**3*x[1]+x[3]*x[2]*q**2*x[1]-q**2*x[2]-q
**2*x[1]**2+x[2]*x[1]*q-x[3]*q*x[2]-x[3]*x[1]**2*q+x[1]*x[3])*q**4*x[1]*x[2]*x[
3]/(-q*x[2]+x[3])/(-1+q*x[1])/(-1+q*x[3])/(x[3]-q)/(-x[1]+x[3])/(-x[1]+q)/(q*x[
1]-x[2]), [0, 0, 0], [q, q*x[1], x[3]]], [-(-1+q*x[1])*(-1+q*x[3])*q**4*x[2]**2
*x[1]*x[3]/(q*x[3]-x[2])/(-1+x[2])**2/(-x[1]+q)/(q*x[1]-x[2])/(x[3]-q), [0, 0, 
0], [q, x[2], q]], [x[3]*x[2]*x[1]*q**4*(q*x[1]-x[2])*(-1+q*x[3])/(-1+q*x[1])/(
x[3]-q)/(q*x[3]-x[2])/(-1+x[2])**2/(-x[1]+q*x[2]), [0, 0, 0], [q*x[2], x[2], q]
], [-x[3]*x[2]*x[1]*q**4*(q-1)*(q+1)*(-x[2]*x[1]+x[3])/(-1+x[2])/(-1+q*x[3])/(x
[3]-q)/(-x[1]+x[3])/(-x[1]+q)/(q*x[3]-x[2]), [0, 0, 0], [q, q*x[3], x[3]]], [-(
q*x[1]-x[2])*(q*x[3]-x[2])*q**4*x[1]*x[3]/(-1+q*x[3])/(-1+q*x[1])/(-x[1]+q*x[2]
)/(-1+x[2])**2/(-q*x[2]+x[3]), [0, 0, 0], [q, 1, q*x[2]]], [x[3]*x[2]*x[1]*q**4
*(-1+q*x[1])*(q*x[3]-x[2])/(q*x[1]-x[2])/(-1+q*x[3])/(-q*x[2]+x[3])/(-1+x[2])**
2/(-x[1]+q), [0, 0, 0], [q, x[2], q*x[2]]], [-(q-1)*(q+1)*(-q**2*x[1]+x[3])*q**
5*x[1]*x[3]*x[2]/(-1+q*x[3])/(-1+q*x[1])/(-x[1]+x[3])/(-q*x[2]+x[3])/(x[3]-q)**
2, [0, 0, 0], [x[3], 1, q]], [-(-q**2*x[1]+x[3])*q**5*x[1]*x[2]/(-1+q*x[3])/(-x
[1]+x[3])/(-q*x[2]+x[3])/(x[3]-q)/(-1+q*x[1]), [0, 0, 0], [x[3], q*x[3], q]], [
-x[2]*x[1]*x[3]**2*q**4*(q-1)**2*(q+1)**2*(-x[1]+q*x[2])/(-x[1]+x[3])/(-q*x[2]+
x[3])/(x[3]-q)/(-x[1]+q)/(-1+q*x[1])/(q*x[1]-x[2])/(-1+q*x[3]), [0, 0, 0], [x[1
], q*x[1], x[3]]], [-x[3]*x[1]*q**4*(q**2*x[1]+q**2*x[2]**2*x[1]-q*x[1]**2-x[2]
*x[1]**2*q-q*x[2]-q*x[2]**2+2*x[2]*x[1])*(q*x[3]-x[2])/(q*x[1]-x[2])/(-1+q*x[3]
)/(-x[1]+q)/(-1+x[2])/(-x[1]+q*x[2])/(-q*x[2]+x[3])/(-1+q*x[1]), [0, 0, 0], [q
, q*x[1], q*x[2]]], [q**4*x[1]*x[2]*x[3]/(-x[1]+q)/(-1+x[2])/(x[3]-q), [0, 0, 1
], [q, 1, q]], [-(q-1)**2*(q+1)**2*(-1+q*x[1])*q**4*x[1]*x[3]**2*x[2]**2/(x[3]-
q)/(q*x[3]-x[2])/(-1+x[2])/(-x[1]+q)/(q*x[1]-x[2])/(-q*x[2]+x[3])/(-1+q*x[3]), 
[0, 0, 0], [q, x[2], x[3]]], [-(q-1)*(q+1)*(-q**3*x[3]+q**2*x[3]**2-x[3]*x[2]*q
**2*x[1]+q**2*x[2]-x[3]*q*x[2]+x[2]*x[1]*q+x[3]**2*x[1]*q-x[1]*x[3])*q**4*x[1]*
x[2]*x[3]/(-1+q*x[3])/(-x[1]+q*x[2])/(-1+q*x[1])/(x[3]-q)/(-x[1]+x[3])/(-x[1]+q
)/(q*x[3]-x[2]), [0, 0, 0], [x[1], q*x[3], q]], [(q-1)*(q+1)*(q**2*x[3]-x[1])*q
**5*x[1]*x[3]*x[2]/(-1+q*x[3])/(-1+q*x[1])/(-x[1]+q)**2/(-x[1]+x[3])/(-x[1]+q*x
[2]), [0, 0, 0], [x[1], 1, q]], [-x[3]*x[2]*x[1]*q**4*(q-1)*(q+1)*(x[3]*x[2]-x[
1])/(-1+x[2])/(-x[1]+q)/(x[3]-q)/(-x[1]+x[3])/(q*x[1]-x[2])/(-1+q*x[1]), [0, 0
, 0], [x[1], q*x[1], q]], [-(q-1)*(q+1)*(-q**2*x[1]+x[3])*q**5*x[1]*x[3]*x[2]/(
-1+q*x[3])/(-1+q*x[1])/(-x[1]+x[3])/(-q*x[2]+x[3])/(x[3]-q)**2, [0, 0, 0], [q, 
1, x[3]]], [-(q*x[1]-x[2])*(q*x[3]-x[2])*q**4*x[1]*x[3]/(-1+q*x[3])/(-1+q*x[1])
/(-x[1]+q*x[2])/(-1+x[2])**2/(-q*x[2]+x[3]), [0, 0, 0], [q*x[2], x[2], q*x[2]]]
, [-x[3]*x[1]**2*x[2]*q**4*(q-1)**2*(q+1)**2*(q*x[3]-x[2])/(-x[1]+q)/(-1+q*x[3]
)/(-x[1]+q*x[2])/(-1+x[2])/(-1+q*x[1])/(q*x[1]-x[2])/(-q*x[2]+x[3]), [0, 0, 0]
, [x[1], x[2], q*x[2]]], [-(q**2*x[3]-x[1])*q**5*x[2]*x[3]/(-1+q*x[3])/(-x[1]+q
)/(-x[1]+x[3])/(-x[1]+q*x[2])/(-1+q*x[1]), [0, 0, 0], [q, q*x[1], x[1]]], [(q-1
)*(q+1)*(q*x[1]-x[2])*q**4*x[1]*x[2]*x[3]/(-1+q*x[1])/(-1+q*x[3])/(x[3]-q)/(q*x
[3]-x[2])/(-1+x[2])/(-x[1]+q*x[2]), [0, 0, 0], [q*x[2], q*x[3], x[3]]], [x[2]*x
[1]*x[3]**2*q**4*(q-1)**2*(q+1)**2*(q*x[1]-x[2])/(q*x[3]-x[2])/(-1+x[2])/(x[3]-
q)/(-1+q*x[1])/(-x[1]+q*x[2])/(-q*x[2]+x[3])/(-1+q*x[3]), [0, 0, 0], [q*x[2], x
[2], x[3]]], [x[2]**2*x[1]**2*x[3]**2*q**4*(q-1)**4*(q+1)**4/(-1+q*x[3])/(-1+q*
x[1])/(-x[1]+q*x[2])/(-q*x[2]+x[3])/(x[3]-q)/(q*x[1]-x[2])/(-x[1]+q)/(q*x[3]-x[
2]), [0, 0, 0], [x[1], x[2], x[3]]], [q**4*x[1]*x[2]*x[3]/(-x[1]+q)/(-1+x[2])/(
x[3]-q), [1, 0, 0], [q, 1, q]], [(-x[1]*x[3]*q*x[2]-q*x[1]*x[3]+2*q*x[3]**2+2*x
[2]*x[1]-x[3]*x[2]-x[3])*q**4*x[1]*x[2]*x[3]/(q*x[3]-x[2])/(-1+q*x[3])/(x[3]-q)
/(-1+x[2])/(-x[1]+q)/(-x[1]+x[3]), [0, 0, 0], [q, q*x[3], q]], [-x[3]*x[1]*q**4
*(-q**2*x[3]-q**2*x[2]**2*x[3]+x[3]**2*q*x[2]+q*x[2]+q*x[2]**2+q*x[3]**2-2*x[3]
*x[2])*(q*x[1]-x[2])/(-1+q*x[1])/(q*x[3]-x[2])/(-1+q*x[3])/(-x[1]+q*x[2])/(x[3]
-q)/(-q*x[2]+x[3])/(-1+x[2]), [0, 0, 0], [q*x[2], q*x[3], q]], [(q-1)*(q+1)*(q*
x[3]-x[2])*q**4*x[1]*x[2]*x[3]/(-x[1]+q)/(-1+q*x[1])/(-1+x[2])/(q*x[1]-x[2])/(-
q*x[2]+x[3])/(-1+q*x[3]), [0, 0, 0], [x[1], q*x[1], q*x[2]]], [-(q*x[1]-x[2])*(
q*x[3]-x[2])*q**4*x[1]*x[3]/(-1+q*x[3])/(-1+q*x[1])/(-x[1]+q*x[2])/(-1+x[2])**2
/(-q*x[2]+x[3]), [0, 0, 0], [q*x[2], 1, q]], [(x[1]*x[3]*q*x[2]+q*x[1]*x[3]-2*x
[3]*x[2]-2*q*x[1]**2+x[2]*x[1]+x[1])*q**4*x[1]*x[2]*x[3]/(q*x[1]-x[2])/(-x[1]+q
)/(-x[1]+x[3])/(-1+x[2])/(x[3]-q)/(-1+q*x[1]), [0, 0, 0], [q, q*x[1], q]], [x[3
]*x[1]**2*x[2]*q**4*(q-1)**2*(q+1)**2*(-q*x[2]+x[3])/(x[3]-q)/(-1+q*x[3])/(-x[1
]+x[3])/(-1+q*x[1])/(-x[1]+q*x[2])/(-x[1]+q)/(q*x[3]-x[2]), [0, 0, 0], [x[1], q
*x[3], x[3]]], [(q-1)**2*(q+1)**2*(-1+q*x[3])*q**4*x[2]**2*x[1]**2*x[3]/(-1+q*x
[1])/(q*x[3]-x[2])/(q*x[1]-x[2])/(x[3]-q)/(-1+x[2])/(-x[1]+q)/(-x[1]+q*x[2]), [
0, 0, 0], [x[1], x[2], q]], [(q-1)*(q+1)*(q**2*x[3]-x[1])*q**5*x[1]*x[3]*x[2]/(
-1+q*x[3])/(-1+q*x[1])/(-x[1]+q)**2/(-x[1]+x[3])/(-x[1]+q*x[2]), [0, 0, 0], [q
, 1, x[1]]], [-2*(-q*x[3]-q*x[1]+2*q**2-q**2*x[2]-3*q**3*x[1]-q**2*x[2]*x[1]**2
+4*q**2*x[3]*x[1]-2*x[3]*x[2]*q**2*x[1]+2*q**2*x[1]**2-3*q**4*x[1]*x[2]*x[3]+2*
q**3*x[2]*x[1]-x[3]**2*x[1]*q-x[3]**2*x[2]*q**2+x[3]*x[2]*x[1]-x[3]*x[1]**2*q+2
*x[3]**2*x[2]*q**3*x[1]-3*q**3*x[3]**2*x[1]+2*x[3]**2*x[1]**2*q**2+2*q**3*x[2]*
x[3]+2*q**2*x[3]**2-3*q**3*x[3]+2*q**3*x[2]*x[1]**2*x[3]-x[3]**2*x[2]*x[1]**2*q
**2-3*q**3*x[1]**2*x[3]+4*q**4*x[1]*x[3])*q**4*x[1]*x[2]*x[3]/(-1+q*x[3])/(-1+q
*x[1])/(-x[1]+q)**2/(-1+x[2])**2/(x[3]-q)**2, [0, 0, 0], [q, 1, q]]}, {}], [{}
, {[1, [0], [q]]}, {}, {[1, [0, 0, 0], [q, 1, q]]}, {}]], [1, 0, 0, 0, 0], [[]
, [x[1]], [x[1], x[2], x[3]], [x[1], x[2], x[3]], []]]
:


#YafeUm(Umb1,F): Given an umbral operator, Umb1, in
#Maple notation, converts it to human notation
#using the symbol F for the function and D for differntiation
YafeUm:=proc(Umb1,F,D)
local T,gu,i,mu,lu,gu1,lu1:

gu:={}:

for i from 1 to nops(Umb1) do
mu:=op(i,Umb1):
gu:=gu union {[mu[2],mu[3]]}:
od:

for i from 1 to nops(gu) do
T[gu[i]]:=0:
od:


for i from 1 to nops(Umb1) do
mu:=op(i,Umb1):
T[[mu[2],mu[3]]]:=normal(T[ [mu[2],mu[3]] ]+mu[1]):
od:

lu:=0:

for i from 1 to nops(gu) do
gu1:=gu[i]:
lu1:=F(op(gu1[2])):

if convert(gu1[1],`+`)<>0 then
 lu1:=lu1*D[op(gu1[1])]:
fi:
lu1:=lu1*factor(T[gu1]):
lu:=lu+lu1:
od:

lu:

end:


#YafeUmSc(UmbSc,F,D): Given an umbral 
#Scheme UmbSc, in
#Maple notation, converts it to human notation
#using the symbols F[1], F[2] for the functions and D for differntiation
#of respective variables
YafeUmSc:=proc(UmbSc,F,D)
local gu1,gu2,gu3,gu4,eq,eq1,n,i,j:
gu1:=UmbSc[1]:
gu2:=UmbSc[2]:
gu3:=UmbSc[3]:
gu4:=UmbSc[4]:
eq:=[]:
n:=nops(gu2):
for i from 1 to n do
eq1:=gu3[i]:

for j from 1 to n do
 eq1:=eq1+YafeUm(gu2[i][j],F[j],D):
od:
eq:=[op(eq),F[i](op(gu4[i]))=eq1]:
od:
eq,gu1:
end:


######################Generalized Stuff



#SAPUmScg(n,x,q,t,w): 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
#followed by the initial values of the corresponding function
#followed by the set of final letters (or states)
#for self-avoiding polygons with the property
#that each vertical cross section has at most n LR-pairs
#here the weight is q^(perimeter)*t^(area)*w^(width)
SAPUmScg:=proc(n,x,q,t,w)
local mu,L,R,A,Z,HATKH,SOF,i1,S,T,U,alist,xlist,gu,gu1,i,j,lu,tu,m,mu1,
yu,zu1:
mu:=Alphabet(n,L,R) union {HATKH,SOF}:
gu:=[]:

for m from 1 to n do

mu1:=Alphabet1(m,L,R):
for i from 1 to nops(mu1) do
xlist:=[seq(x[i1],i1=1..2*m-1)]:
alist:=[seq(A[i1],i1=1..2*m-1)]:
S[mu1[i]]:=xlist:
 U[mu1[i]]:=0:

lu:=PreUmbrag(mu1[i],L,R,n,x,q,t,w,A,Z,HATKH,SOF):

for j from 1 to nops(mu) do
T[mu1[i],mu[j]]:=ToUmbra(coeff(lu,Z[mu[j]],1),xlist,alist):
od:


od:

od:

mu1:={HATKH,SOF}:
for i from 1 to nops(mu1) do
xlist:=[]:
alist:=[]:
S[mu1[i]]:=xlist:

if mu1[i]=HATKH then
 U[mu1[i]]:=1:
else
 U[mu1[i]]:=0:
fi:

lu:=PreUmbrag(mu1[i],L,R,n,x,q,t,w,A,Z,HATKH,SOF):

for j from 1 to nops(mu) do
T[mu1[i],mu[j]]:=ToUmbra(coeff(lu,Z[mu[j]],1),xlist,alist):
od:
od:


gu:=[]:
tu:=[]:

mu:=convert(mu,list):

yu:=[]:

for i from 1 to nops(mu) do
 tu:=[op(tu),S[mu[i]]]:
 yu:=[op(yu),U[mu[i]]]:
 if mu[i]=SOF then
  zu1:={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:


 

#PreUmbrag(LETTER,L,R,n,x,q,t,w,A,Z,HATKH,SOF)
#Given a letter LETTER, in terms of legal
#parantheses using L and R
#and an integer n specifying the max. number of
#LR-pairs, and a letter Z
#and a symbol A, finds the
#pre-umbra (i.e. the image of a generic monomial
#x[1]^A[1]*x[2]^A[2]*...*x[r]^A[r] where r is nops(LETTER)-1
#as a linear combination of Z[LETTER1], where LETTER1
#ranges over all followers of LETTER
#the weight is q^perimter*t^area*w^(width)

PreUmbrag:=proc(LETTER,L,R,n,x,q,t,w,A,Z,HATKH,SOF)
local gu,mu,C,i,n1,lu,tu,i1:

if LETTER=SOF then
RETURN(0):
fi:

if LETTER=HATKH then
 gu:=0:
 for i from 1 to n do
  lu:=[seq(op([L,R]),i1=1..i)]:
  tu:=q^3*t*x[1]/(1-q*t*x[1]):
   for i1 from 1 to i-1 do
    tu:=tu*x[2*i1]/(1-x[2*i1])*q^3*x[2*i1+1]/(1-q*x[2*i1+1]):
   od:
   gu:=gu+tu*w*Z[lu]:
 od:
RETURN(gu):
fi:
  
mu:=Followers(LETTER,L,R,C,n):

gu:=0:

for i from 1 to nops(mu) do
gu:=gu+w*APUg(LETTER,mu[i],L,R,C,x,q,t,A)*Z[OutLetter(mu[i],L,R,C)]:
od:

n1:=nops(LETTER)/2:
if LETTER=[L,seq(op([L,R]),i=1..n1-1),R] then
lu:=0:
 for i from 1 to n1 do
  lu:=lu+A[2*i-1]:
 od:

  gu:=gu+Z[SOF]*q^lu:
fi:

expand(gu):
end:


#APUg(InpLet,OutPreLet,L,R,C,x,q,t,A): Given an input letter from the left
#InpLet, and an output Pre-letter from the right,
#OutPreLet, finds the Atomic-Pre-Umbra
#with respect to the "length-variables" x[1],x[2],...
#and q is the variable that keeps track of the length
#of the partial SAP
#the weight is q^(perimeter)*t^(area)
APUg:=proc(InpLet,OutPreLet,L,R,C,x,q,t,A)
local w,i,gu1,gu2,gu3,gu4,mu,gu,n,T,i1,lu,T0,Tinf,xinf,lub,lut:
n:=nops(InpLet):
gu1:=OutPreLet[1]:
gu2:=OutPreLet[2]:
gu3:=OutPreLet[3]:
gu4:=OutPreLet[4]:
w:=[]:
for i from 1 to nops(gu1) do
if gu1[i]=C then
 w:=[op(w),C]:
else
 w:=[op(w),gu1[i][1]]:
fi:
od:

for i from 1 to nops(gu2) do
 T[gu2[i][1]]:=0:
od:

for i from 1 to nops(gu3) do
 T[gu3[i][1]]:=gu4[i]:
od:

gu:=[[x[0],seq(op([q*x[2*i1-1],x[2*i1]]),i1=1..T[0])]]:
mu:=2*T[0]:
for i from 1 to n do
if w[i]=C then
 gu:=[op(gu),[x[mu],seq(op([q*x[mu+2*i1+1],x[mu+2*i1+2]]),i1=0..T[i]-1)]]:
 mu:=mu+2*T[i]:
else
 gu:=[op(gu),[x[mu+1],seq(op([q*x[mu+2*i1+2],x[mu+2*i1+3]]),i1=0..T[i]-1)]]:
 mu:=mu+2*T[i]+1:
fi:
od:
gu:


for i from 1 to nops(gu) do
T[i-1]:=gu[i]:
od:

for i from 1 to nops(gu1) do
if gu1[i][1]<>C then
 if gu1[i][2]=-1 then
 T[i-1]:=[op(T[i-1]),q*T[i][1]]:
 elif gu1[i][2]=1 then
 T[i]:=[q*T[i-1][nops(T[i-1])],op(T[i])]:
 fi:
fi:
od:

gu:=[seq(T[i],i=0..nops(gu1))]:

lu:=q^(nops(OutLetter(OutPreLet,L,R,C))):

for i from 1 to nops(gu2) do
lu:=lu*q^(A[gu2[i][1]]):
od:

T0:=[op(2..nops(T[0]),T[0])]:
Tinf:=[op(1..nops(T[nops(gu1)])-1,T[nops(gu1)])]:
xinf:=op(nops(T[nops(gu1)]),T[nops(gu1)]):

for i1 from 1 to nops(T0) do
lu:=lu*T0[i1]/(1-T0[i1]):
od:

for i from 1 to nops(gu1)-1 do
 lu:=lu*PzA(T[i],A[i]):
od:

for i1 from 1 to nops(Tinf) do
lu:=lu*Tinf[i1]/(1-Tinf[i1]):
od:

lu:=normal(lu):
lu:=normal(subs({x[0]=1},lu)):
if expand(subs(xinf=1,denom(lu)))=0 then
  if expand(subs(xinf=1,numer(lu)))<>0 then
   ERROR(`Something is wrong`):
  fi:

 lut:= subs(xinf=1,diff(numer(lu),xinf)):
 lub:= subs(xinf=1,diff(denom(lu),xinf)):
 lu:=normal(lut/lub):
else
lu:=normal(subs({xinf=1},lu)):
fi:

for i from 1 by 2 to n do
lu:=subs(x[i]=t*x[i],lu):
od:
lu:

end:


# ApplyUmScg(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
ApplyUmScg:=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(normal(lu1[i1]),q,i))]:
od:
lu:=gu:
print(`lu is`,lu):
fu:=[]:
for i1 from 1 to nops(lu) do
fu:=[op(fu),coeff(expand(lu[i1]),q,i-1)]:
od:
print(`fu is`,fu):
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))}:
print(`mu before is`, normal(mu)):
mu:=subs(ku,normal(mu)):
print(`ku is`, ku):
print(`mu after is`, normal(mu)):

Sid:=[op(Sid),normal(mu)]:

lprint(Sid):
od:

lu,Sid:

end:


#########Stuff for SAPS with restricted extensions

#FollowersR(LET,L,R,C,n,Rset): All the ways to follow
#the letter LET if the maximal number of pairs  of
#parantheses is n, where the extensions should be
#in the set Rset
FollowersR:=proc(LET,L,R,C,n,Rset)
local gu,i,mu,j,w:

mu:=Followers(LET,L,R,C,n):

gu:={}:

for i from 1 to nops(mu) do
w:=mu[i][1]:

if member([seq(w[j][2],j=1..nops(w))],Rset) then
 gu:=gu union {mu[i]}:
fi:

od:
gu:
end:


#PreUmbragR(LETTER,L,R,n,x,q,t,w,A,Z,HATKH,SOF,Rset)
#Like PreUmbrag(LETTER,L,R,n,x,q,t,w,A,Z,HATKH,SOF)
#but for the resticted set of SAPS such that the continuations
#belong to Rset
PreUmbragR:=proc(LETTER,L,R,n,x,q,t,w,A,Z,HATKH,SOF,Rset)
local gu,mu,C,i,n1,lu,tu,i1:

if LETTER=SOF then
RETURN(0):
fi:

if LETTER=HATKH then
 gu:=0:
 for i from 1 to n do
  lu:=[seq(op([L,R]),i1=1..i)]:
  tu:=q^3*t*x[1]/(1-q*t*x[1]):
   for i1 from 1 to i-1 do
    tu:=tu*x[2*i1]/(1-x[2*i1])*q^3*x[2*i1+1]/(1-q*x[2*i1+1]):
   od:
   gu:=gu+tu*w*Z[lu]:
 od:
RETURN(gu):
fi:
  
mu:=FollowersR(LETTER,L,R,C,n,Rset):

gu:=0:

for i from 1 to nops(mu) do
gu:=gu+w*APUg(LETTER,mu[i],L,R,C,x,q,t,A)*Z[OutLetter(mu[i],L,R,C)]:
od:

n1:=nops(LETTER)/2:
if LETTER=[L,seq(op([L,R]),i=1..n1-1),R] then
lu:=0:
 for i from 1 to n1 do
  lu:=lu+A[2*i-1]:
 od:

  gu:=gu+Z[SOF]*q^lu:
fi:

expand(gu):
end:



#SAPUmScgR(n,x,q,t,w,Rset): 
#Like SAPUmScgR(n,x,q,t,w), but only for SAPS for which the
#continuations belong to the set Rset
SAPUmScgR:=proc(n,x,q,t,w,Rset)
local mu,L,R,A,Z,HATKH,SOF,i1,S,T,U,alist,xlist,gu,gu1,i,j,lu,tu,m,mu1,
yu,zu1:
mu:=Alphabet(n,L,R) union {HATKH,SOF}:
gu:=[]:

for m from 1 to n do

mu1:=Alphabet1(m,L,R):

for i from 1 to nops(mu1) do
xlist:=[seq(x[i1],i1=1..2*m-1)]:
alist:=[seq(A[i1],i1=1..2*m-1)]:
S[mu1[i]]:=xlist:
 U[mu1[i]]:=0:

lu:=PreUmbragR(mu1[i],L,R,n,x,q,t,w,A,Z,HATKH,SOF,Rset):

for j from 1 to nops(mu) do
T[mu1[i],mu[j]]:=ToUmbra(coeff(lu,Z[mu[j]],1),xlist,alist):
od:


od:

od:

mu1:={HATKH,SOF}:
for i from 1 to nops(mu1) do
xlist:=[]:
alist:=[]:
S[mu1[i]]:=xlist:

if mu1[i]=HATKH then
 U[mu1[i]]:=1:
else
 U[mu1[i]]:=0:
fi:

lu:=PreUmbrag(mu1[i],L,R,n,x,q,t,w,A,Z,HATKH,SOF):

for j from 1 to nops(mu) do
T[mu1[i],mu[j]]:=ToUmbra(coeff(lu,Z[mu[j]],1),xlist,alist):
od:
od:


gu:=[]:
tu:=[]:

mu:=convert(mu,list):

yu:=[]:

for i from 1 to nops(mu) do
 tu:=[op(tu),S[mu[i]]]:
 yu:=[op(yu),U[mu[i]]]:
 if mu[i]=SOF then
  zu1:={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:


#SAP1UmScR(x,q,Rset): 
#A simplified Umbral Scheme in the form
#[INI,Rota_Operator]
#for SAPS with only two horizontal edges per vertical 
#cross-section, in which the follow-up rules
#are according to Rset, where 1 means up, 0 means level,
#and -1 means down
#For example SAP1UmScR(x,q,{[0,0],[0,1],[1,0],[1,1]}) 
#gives the Umbral Scheme for suchs SAPS in which
#you never go down from a horizontal edge, but have
#a weakly upwards trend
SAP1UmScR:=proc(x,q,Rset)  
local gu,gu2,gu3,gu4,INI,j,i,GU,GU1,lu: 
gu:=SAPUmScgR(1,x,q,1,1,Rset):
gu2:=gu[2]:
gu3:=gu[3]:
gu4:=gu[4]:

for j from 1 to nops(gu3) while gu3[j]=0 do od:
for i from 1 to nops(gu4) while gu4[i]<>[x[1]] do od:
GU:=gu2[i][i]:
INI:=gu2[i][j][1][1]:
GU:=subs(x[1]=x,GU):
INI:=subs(x[1]=x,INI):
GU1:={}:
for i from 1 to nops(GU) do
lu:=GU[i]:
GU1:=GU1 union {[lu[1],lu[2][1],lu[3][1]]}:
od:
[INI,GU1,x]:
end:


#FunEq(SUE,A): Given a Simplified Umbral Scheme
#(with one function), and a symbol A, finds
#the Functional Equation in a Human Readable Form
FunEq:=proc(SUE,A)
local eq,INI,GU,x,i,lu:
INI:=SUE[1]: GU:=SUE[2]: x:=SUE[3]:
eq:=INI:

for i from 1 to nops(GU) do
lu:=GU[i]:
eq:=eq-lu[1]*subs(x=lu[3],xdiff1(A(x),x,lu[2])):
od:
eq:
end:

#Maylis(x,q,Rset,A): 
#The Functional Equation, in Human Form, for the
#generating function, let's call it A(x) for SAPS
#with exactly two horizontal edges in every vertical
#cross-section, with the interfaces restriced by the
#set of pairs Rset
#according to the weight q^parameter x^(width of last vertical wall)
Maylis:=proc(x,q,Rset,A):
FunEq(SAP1UmScR(x,q,Rset),A)  : 
end:



#SSUE(SUE,q): Given a Simplified Umbral Equation, SUE, and
# a variable q, that is a linear combinatrion of
#A(x),A(q), and qA'(q), solves for A(q)
SSUE:=proc(SUE,q)
local x,eq,var,gu,INI,GU,i,lu,Ax,Aq,Aq1,gu1,Ax1,gu2,mu,yeter,
pita,pita1,ku,mu1:
INI:=SUE[1]: GU:=SUE[2]: x:=SUE[3]:

gu:=Ax-INI:
for i from 1 to nops(GU) do
lu:=GU[i]:
if lu[3]=x then
 if lu[2]=0 then
  gu:=gu-lu[1]*Ax:
  else
   print(`This Simplified Umbral Equation Is not treated by this procedure`):
   RETURN(0):
 fi:

elif lu[3]=q then
 if lu[2]=0 then
  gu:=gu-lu[1]*Aq:
 elif lu[2]=1 then
  gu:=gu-lu[1]*Aq1:
   else
  print(`This Simplified Umbral Equation Is not treated by this procedure`):
  RETURN(0):
 fi:
else
  print(`This Simplified Umbral Equation Is not treated by this procedure`):
  RETURN(0):
fi:
od:

gu:=expand(numer(normal(gu))):

gu1:=subs(x=q,gu):
gu1:=expand(subs(Ax=Aq,gu1)):
if coeff(gu1,Aq,1)<>0 then
 if coeff(gu,Aq1,1)=0 then
 RETURN(solve(gu1,Aq)):
 else
 gu2:=x*diff(coeff(gu,Ax,1),x)*Ax+coeff(gu,Ax,1)*Ax1+
x*diff(expand(gu-coeff(gu,Ax,1)*Ax),x):
gu2:=expand(gu2):
gu2:=subs(x=q,gu2):
gu2:=subs({Ax=Aq,Ax1=Aq1},gu2):
eq:={gu1,gu2}:
var:={Aq,Aq1}:
var:=solve(eq,var):

if var=NULL then
print(`No solution found`):
RETURN(0):
fi:

RETURN(normal(subs(var,Aq))):
fi:
fi:

if coeff(gu1,Aq,1)=0 then

mu:=coeff(gu,Ax,1):
mu1:=factor(numer(factor(expand(mu/coeff(mu,x,degree(mu,x)))))):
yeter:=expand(gu-mu*Ax):

if degree(mu1,x)>2 and type(factor(mu1),`+`) then
print(`I rather not solve irreducible equations of degreee >=3 right now`):
RETURN(0):
fi:
pita:={solve(mu,x)}:
pita1:={}:
for i from 1 to nops(pita) do
if (type(numer(normal(pita[i])),polynom) and type(denom(normal(pita[i])),
polynom)) or
(expand(subs(q=0,numer(normal(pita[i]))))=0 and 
expand(subs(q=0,denom(normal(pita[i]))))=0)
 then
 pita1:=pita1 union {pita[i]}:
fi:
od:

var:={Aq,Aq1}:
eq:={}:
for i from 1 to nops(pita1) do
ku:=expand(numer(normal(subs(x=pita1[i],yeter)))):
if (coeff(ku,Aq)<>0 or coeff(ku,Aq1)<>0) then
eq:=eq union {ku}:
fi:
 
od:


var:=solve(eq,var):

if var=NULL then
print(`No solution found`):
RETURN(0):
fi:

RETURN(normal(subs(var,Aq))):

else
print(`I can't solve it`):
RETURN(0):
fi:

end:



#SSUEx(SUE,q,x): Given a Simplified Umbral Equation, SUE, and
# a variable q, that is a linear combinatrion of
#A(x),A(q), and qA'(q), solves for A(q)
#it also returns A(x)
SSUEx:=proc(SUE,q,x)
local eq,var,gu,INI,GU,i,lu,Ax,Aq,Aq1,gu1,Ax1,gu2,mu,yeter,pita,
pita1,AqS,ku:
INI:=SUE[1]: GU:=SUE[2]: x:=SUE[3]:

gu:=Ax-INI:
for i from 1 to nops(GU) do
lu:=GU[i]:
if lu[3]=x then
 if lu[2]=0 then
  gu:=gu-lu[1]*Ax:
  else
   print(`This Simplified Umbral Equation Is not treated by this procedure`):
  RETURN(0):
 fi:

elif lu[3]=q then
 if lu[2]=0 then
  gu:=gu-lu[1]*Aq:
 elif lu[2]=1 then
  gu:=gu-lu[1]*Aq1:
   else
  print(`This Simplified Umbral Equation Is not treated by this procedure`):
  RETURN(0):
 fi:
else
  print(`This Simplified Umbral Equation Is not treated by this procedure`):
  RETURN(0):
fi:
od:
gu:=expand(numer(normal(gu))):
gu1:=subs(x=q,gu):
gu1:=expand(subs(Ax=Aq,gu1)):

if coeff(gu1,Aq,1)<>0 then
 if coeff(gu,Aq1,1)=0 then
  AqS:=normal(solve(gu1,Aq)):
  RETURN(AqS, normal(solve( subs(Aq=AqS,gu),Ax))):
 else
 gu2:=x*diff(coeff(gu,Ax,1),x)*Ax+coeff(gu,Ax,1)*Ax1+
x*diff(expand(gu-coeff(gu,Ax,1)*Ax),x):
gu2:=expand(gu2):
gu2:=subs(x=q,gu2):
gu2:=subs({Ax=Aq,Ax1=Aq1},gu2):
eq:={gu1,gu2}:
var:={Aq,Aq1}:
var:=solve(eq,var):

if var=NULL then
print(`No solution found`):
RETURN(0):
fi:

RETURN(normal(subs(var,Aq)),normal(solve(subs(var,gu),Ax))):
fi:
fi:

if coeff(gu1,Aq,1)=0 then

mu:=coeff(gu,Ax,1):
yeter:=expand(gu-mu*Ax):
pita:={solve(mu,x)}:
if pita={} then
ERROR(`pita is empty`):
fi:

pita1:={}:
for i from 1 to nops(pita) do

 if (type(numer(pita[i]),polynom) and type(denom(pita[i]),polynom)) or
(subs(q=0,numer(pita[i]))=0 and subs(q=0,denom(pita[i]))=0 ) then
 pita1:=pita1 union {pita[i]}:
fi:
od:

var:={Aq,Aq1}:
eq:={}:

for i from 1 to nops(pita1) do
ku:=expand(numer(normal(subs(x=pita1[i],yeter)))):
if (coeff(ku,Aq)<>0 or coeff(ku,Aq1)<>0) then
eq:=eq union {ku}:
fi:

od:
var:=solve(eq,var):
if var=NULL then
ERROR(`Can't solve it`):
fi:
RETURN(normal(subs(var,Aq)),normal(solve(subs(var,gu),Ax))):

else
print(`I can't solve it`):
RETURN(0):
fi:

end:


#Xavier(Rset,q): The Generating Function for column-convex
#polyominoes (equivalenty, SAPs with at most two horizontal
#edges per vertical cross-section) according to parameter,
#where interfaces are constrained to belong to the set
#Rset, for example, try: Xavier({[1,1]},q);
Xavier:=proc(Rset,q) local x: 
if Rset={} then RETURN(1): fi:
SSUE(SAP1UmScR(x,q,Rset),q): end:

#Xavierx(Rset,q,x): Like Xavier(Rset,q) but also
#returns the more detailed generating function
#with wt=q^perimerter*x^(length of rightmost wall)
Xavierx:=proc(Rset,q,x): if Rset={} then RETURN(1): fi:
SSUEx(SAP1UmScR(x,q,Rset),q,x): end:


#Delest(RsetL,q): The Generating Function for column-convex
#polyominoes (equivalenty, SAPs with at most two horizontal
#edges per vertical cross-section) according to parameter,
#where it firsts goes according to the
#interfaces in RsetL[1], then it either quits or starts
#going according to the
#interfaces in RsetL[2], then it quits, or starts going
#according to the interfaces in RsetL[3], and so on
#It is used in ProveDelestViennot
Delest:=proc(RsetL,q)
local lu,TOT,GU,INI,mu,x,LU,INI1,i,i1:

GU:=SAP1UmScR(x,q,RsetL[1]):
mu:=SSUEx(GU,q,x):
TOT:=mu[1]:
INI:=mu[2]:

for i from 2 to nops(RsetL) do
lu:=RsetL[i] minus RsetL[i-1]:
LU:=SAP1UmScR(x,q,lu)[2]:
INI1:=0:
for i1 from 1 to nops(LU) do
INI1:=INI1+LU[i1][1]*subs(x=LU[i1][3],xdiff1(INI,x,LU[i1][2])):
od:
INI1:=normal(INI1):
GU:=[INI1,SAP1UmScR(x,q,RsetL[i])[2],x]:
mu:=SSUEx(GU,q,x):
TOT:=normal(TOT+mu[1]):
INI:=mu[2]:
od:
TOT:
end:


#Delest1(RsetL,q): The Generating Function for column-convex
#polyominoes (equivalenty, SAPs with at most two horizontal
#edges per vertical cross-section) according to parameter,
#where it firsts goes according to the
#interfaces in RsetL[1], then it  starts
#going according to the
#interfaces in RsetL[2], then it starts going
#according to the interfaces in RsetL[3], and so on
#For example, try: 
#Delest1([{[-1,0],[-1,1],[0,0],[0,1]},
#{[0,0],[0,1],[1,0],[1,1]},
#{[0,0],[0,-1],[1,0],[1,-1]}],q);,
Delest1:=proc(RsetL,q)
local lu,GU,INI,mu,x,LU,INI1,i,i1:

GU:=SAP1UmScR(x,q,RsetL[1]):
mu:=SSUEx(GU,q,x):
INI:=mu[2]:

for i from 2 to nops(RsetL) do
lu:=RsetL[i] minus RsetL[i-1]:
LU:=SAP1UmScR(x,q,lu)[2]:
INI1:=0:
for i1 from 1 to nops(LU) do
INI1:=INI1+LU[i1][1]*subs(x=LU[i1][3],xdiff1(INI,x,LU[i1][2])):
od:
INI1:=normal(INI1):
GU:=[INI1,SAP1UmScR(x,q,RsetL[i])[2],x]:
mu:=SSUEx(GU,q,x):
INI:=mu[2]:
od:
mu[1]:
end:

#DV(q): The Delest-Viennot Generating Function for the
#number of convex polygons according to perimeter
DV:=proc(q)
local dv,gu:
dv:=[{[-1,0],[-1,1],[0,0],[0,1]},
{[0,0],[0,1],[1,0],[1,1]},
{[0,0],[0,-1],[1,0],[1,-1]}];
gu[1]:=Delest1([dv[1]],q):
gu[2]:=Delest1([dv[1],dv[2]],q):
gu[3]:=Delest1(dv,q):
gu[4]:=Delest1([dv[1],dv[3]],q):
gu[5]:=Delest1([dv[1],dv[2] intersect dv[3]],q):
gu[6]:=Delest1([dv[1],dv[2] intersect dv[3],dv[3]],q):
normal(gu[1]+2*(gu[2]-gu[5])+2*(gu[3]-gu[6])+gu[4]):
end:


#ProveDelestViennot(): Proves the Delest-Viennot Generating Function
#(equivalent to the explicit formula a(2n+8)=(2n+11)4^n-4(2n+1)!/n!^2,
#for n>=0, and a(4)=1, a(6)=2)
ProveDelestViennot:=proc()
local dv,gu,CONJ,q,n,t0,lu:
t0:=time():
CONJ:=q^4+2*q^6+
sum((2*n+11)*4^n*q^(2*n+8),n=0..infinity)-
sum(4*(2*n+1)*binomial(2*n,n)*q^(2*n+8),n=0..infinity):
dv:=[{[-1,0],[-1,1],[0,0],[0,1]},
{[0,0],[0,1],[1,0],[1,1]},
{[0,0],[0,-1],[1,0],[1,-1]}];
gu[1]:=Delest1([dv[1]],q):
gu[2]:=Delest1([dv[1],dv[2]],q):
gu[3]:=Delest1(dv,q):
gu[4]:=Delest1([dv[1],dv[3]],q):
gu[5]:=Delest1([dv[1],dv[2] intersect dv[3]],q):
gu[6]:=Delest1([dv[1],dv[2] intersect dv[3],dv[3]],q):
print(`The Delest-Viennot celebrated result that`):
print(`the number of convex polyominoes with perimeter 2n+8`):
print(` equals (2*n+11)*4^n-4*(2*n+1)!/n!^2  is indeed`):
lu:=gu[1]+2*(gu[2]-gu[5])+2*(gu[3]-gu[6])+gu[4]:
lu:=normal(lu):
print(evalb(0=normal(CONJ-
(gu[1]+2*(gu[2]-gu[5])+2*(gu[3]-gu[6])+gu[4])))):
print(`The whole proof took`, time()-t0, `seconds of CPU time`):

end:


StrictDV:=proc(q)
local dv,gu:
dv:=[{[-1,1]},
{[1,1]},
{[1,-1]}];
gu[1]:=Delest1([dv[1]],q):
gu[2]:=Delest1([dv[1],dv[2]],q):
gu[3]:=Delest1(dv,q):
gu[4]:=Delest1([dv[1],dv[3]],q):
normal(gu[1]+2*gu[2]+2*gu[3]+gu[4]):
end:



#DoThemAllT(TOT,q): Tries to Find the generating funtion of
#All Consistent classes of Column-Convex polyominoes
#with the further restriction that the interfaces
#obey the set Rset, for all possible (inequivalent)
#subsets of the full gamut (for column-convex polyominoes)
#
DoThemAllT:=proc(TOT,q)
local gu,mu,t0,lu,mu1,i,t00,SUC,FAI:
t00:=time():
with(combinat):
SUC:=0: FAI:=0:
gu:=powerset(TOT) minus {{}}:
while gu<>{} do
mu:=op(1,gu):
t0:=time():
lu:=Xavierterse(mu,q):
if lu<>0 then
SUC:=SUC+1:
print(`The generating function for Column-Convex Polyominoes according`):
print(`to perimeter, where the interfaces are restricted to lie in`):
print(mu):
print(` is to equal `):
print(lu):
print(`It took`, time()-t0, `seconds of CPU time`):
print(`The first 20 terms are `):
lu:=taylor(lu,q=0,45):
print([seq(coeff(lu,q,2*i),i=2..20)]):
else
FAI:=FAI+1:
fi:

mu1:={seq([-mu[i][1],-mu[i][2]],i=1..nops(mu))}:
gu:=gu minus {mu,mu1}:
od:
print(`I succeeded in `, SUC, `cases and failed in`, FAI, ` cases `):
print(`Everything took`, time()-t00, `seconds of CPU time `):
end:

#DoThemAll(q): Tries to Find the generating funtion of
#All Consistent classes of Column-Convex polyominoes
#with the further restriction that the interfaces
#obey the set Rset, for all possible (inequivalent)
#subsets of the full gamut (for column-convex polyominoes)
#
DoThemAll:=proc(q):
DoThemAllT({[-1,-1],[-1,0],[-1,1],[0,-1],[0,0],[0,1],
[1,-1],[1,0],[1,1]},q):
end:


#SSUEterse(SUE,q): Terse version of SSUE
SSUEterse:=proc(SUE,q)
local x,eq,var,gu,INI,GU,i,lu,Ax,Aq,Aq1,gu1,Ax1,gu2,mu,yeter,
pita,pita1,ku,mu1:
INI:=SUE[1]: GU:=SUE[2]: x:=SUE[3]:

gu:=Ax-INI:
for i from 1 to nops(GU) do
lu:=GU[i]:
if lu[3]=x then
 if lu[2]=0 then
  gu:=gu-lu[1]*Ax:
  else
  RETURN(0):
 fi:

elif lu[3]=q then
 if lu[2]=0 then
  gu:=gu-lu[1]*Aq:
 elif lu[2]=1 then
  gu:=gu-lu[1]*Aq1:
   else
  RETURN(0):
 fi:
else
  RETURN(0):
fi:
od:

gu:=expand(numer(normal(gu))):

gu1:=subs(x=q,gu):
gu1:=expand(subs(Ax=Aq,gu1)):
if coeff(gu1,Aq,1)<>0 then
 if coeff(gu,Aq1,1)=0 then
 RETURN(solve(gu1,Aq)):
 else
 gu2:=x*diff(coeff(gu,Ax,1),x)*Ax+coeff(gu,Ax,1)*Ax1+
x*diff(expand(gu-coeff(gu,Ax,1)*Ax),x):
gu2:=expand(gu2):
gu2:=subs(x=q,gu2):
gu2:=subs({Ax=Aq,Ax1=Aq1},gu2):
eq:={gu1,gu2}:
var:={Aq,Aq1}:
var:=solve(eq,var):

if var=NULL then
RETURN(0):
fi:

RETURN(normal(subs(var,Aq))):
fi:
fi:

if coeff(gu1,Aq,1)=0 then

mu:=coeff(gu,Ax,1):
mu1:=factor(numer(factor(expand(mu/coeff(mu,x,degree(mu,x)))))):
yeter:=expand(gu-mu*Ax):

if degree(mu1,x)>2 and type(factor(mu1),`+`) then
RETURN(0):
fi:
pita:={solve(mu,x)}:
pita1:={}:
for i from 1 to nops(pita) do
if (type(numer(normal(pita[i])),polynom) and type(denom(normal(pita[i])),
polynom)) or
(expand(subs(q=0,numer(normal(pita[i]))))=0 and 
expand(subs(q=0,denom(normal(pita[i]))))=0)
 then
 pita1:=pita1 union {pita[i]}:
fi:
od:

var:={Aq,Aq1}:
eq:={}:
for i from 1 to nops(pita1) do
ku:=expand(numer(normal(subs(x=pita1[i],yeter)))):
if (coeff(ku,Aq)<>0 or coeff(ku,Aq1)<>0) then
eq:=eq union {ku}:
fi:
 
od:


var:=solve(eq,var):

if var=NULL then
RETURN(0):
fi:

RETURN(normal(subs(var,Aq))):

else
RETURN(0):
fi:

end:
Xavierterse:=proc(Rset,q) local x: 
if Rset={} then RETURN(1): fi:
SSUEterse(SAP1UmScR(x,q,Rset),q): end:

#########End Stuff for SAPS with restricted extensions