print(`AARON: A small Maple package accompanying an article by`):
print(`Robertson, Saracino, and Zeilberger`):
print(``):
print(`You will need the Maple package WILF written by Zeilberger`):
print(``):
print(`For a list of the procedures and the syntax, type: ez(); `):
read WILF:
ez:=proc(): print(`T(per), ImageT(n), ImageS(n),Cata(n),S(cata) `):
print(`CheckST(n), FixedPts(per), DiagPts, Anx(n,x), Bnx(n,x),anx(n,x)`):
print(`Cn(n), Enrx(n,r,x), enrx(n,r,x), numbrrp(pat,n,k,w) `):
print(`bnk(n,k), bnkchk(n), bnk2(n,k), bnk3(n,k),catchk(n),fnchk(n)`):
end:


#T(per): the mapping T from 321-avoiding permutations of length
#n and Catalan sequences a1<=a2<=...<=an such that ai<=i
T:=proc(per)
local i,n,per1:
n:=nops(per):
if n=0 then
 RETURN([]):
fi:
if n=1 then
RETURN([1]):
fi:

for i from 1 to nops(per) while per[i]<>n do
od:
if per[n]=n-1 then
per1:=redu([op(1..n-1,per)]):
else
per1:=redu([op(1..i-1,per),op(i+1..n,per)]):
fi:
[op(T(per1)),i]:

end:


S:=proc(cata)
local cata1,i,j,n,per1:
n:=nops(cata):
if cata=[] then
RETURN([]):
fi:
if cata=[1] then
 RETURN([1]):
fi:
i:=cata[n]:
j:=cata[n-1]:
cata1:=[op(1..n-1,cata)]:
per1:=S(cata1):
if i=j then
RETURN([op(1..i-1,per1),n,op(i+1..n-1,per1),n-1]):
else
RETURN([op(1..i-1,per1),n,op(i..n-1,per1)]):
fi:
end:



#Cata(n): The Cata sequences of length n
Cata:=proc(n)
local lu,i,gu,lu1,j,akha:
option remember:
if n=0 then
 RETURN({[]}):
fi:
if n=1 then
RETURN({[1]}):
fi:
lu:=Cata(n-1):
gu:={}:
for i from 1 to nops(lu) do
lu1:=lu[i]:
akha:=lu1[nops(lu1)]:
gu:=gu union {seq([op(lu1),j],j=akha..n)}:
od:
gu:
end:

ImageT:=proc(n) local lu,i: 
lu:=Wilf(n,{[3,2,1]}):{seq(T(lu[i]),i=1..nops(lu))}:end:
ImageS:=proc(n) local lu,i: 
lu:=Cata(n):{seq(S(lu[i]),i=1..nops(lu))}:end:


#CheckST(n): check whether ST(per)=per for every 321-avoiding
#permutation of length n
CheckST:=proc(n)
local gu,i:
gu:=Wilf(n,{[3,2,1]}):
for i from 1 to nops(gu) do
if not (S(T(gu[i]))=gu[i]  and FixedPts(gu[i])=DiagPts(T(gu[i])) ) then
 RETURN(false):
fi:
od:
true:
end:


#CheckTS(n): check whether TS(cata)=cata for every Catalan sequence
#of length n
CheckTS:=proc(n)
local gu,i:
gu:=Cata(n):
for i from 1 to nops(gu) do
if not (T(S(gu[i]))=gu[i] and FixedPts(S(gu[i]))=DiagPts(gu[i])) then
 RETURN(false):
fi:
od:
true:
end:


#FixedPts(per): the set of fixed points of a permutation
FixedPts:=proc(per)
local gu,i:
gu:={}:
for i from 1 to nops(per) do
if per[i]=i then
 gu:=gu union {i}:
fi:
od:
gu:
end:

#DiagPts(cata): the set of i such that cata[i]=i and cata[i+1]=i+1
#(i<nops(cata)) and cata[n]=n
DiagPts:=proc(cata)
local i,gu,n:
n:=nops(cata):
gu:={}:
for i from 1 to n-1 do
if cata[i]=i and cata[i+1]=i+1 then
 gu:=gu union {i}:
fi:
od:
if cata[n]=n then
gu:=gu union {n}:
fi:
gu:
end:

#Anx(n,x): The generating function for 321-avoiding permutations
#from first principle
Anx:=proc(n,x)
local gu,i,lu:
lu:=Wilf(n,{[3,2,1]}):
gu:=0:
for  i from 1 to nops(lu) do
gu:=gu+x^nops(FixedPts(lu[i])):
od:
gu:
end:

#Bnx(n,x): The generating function for 132-avoiding permutations
#from first principle
Bnx:=proc(n,x)
local gu,i,lu:
lu:=Wilf(n,{[1,3,2]}):
gu:=0:
for  i from 1 to nops(lu) do
gu:=gu+x^nops(FixedPts(lu[i])):
od:
gu:
end:

#anx(n,x): Anx(n,x) using the generating function
anx:=proc(n,x)
local t,lu,psi:
psi:=(1-sqrt(1-4*t))/(2*t):
lu:=(1-x+psi)/((2-x)+(x-1)^2*t):
lu:=taylor(lu,t=0,n+3):
normal(coeff(lu,t,n)):
end:

Cn:=proc(n):binomial(2*n,n)/(n+1):end:

#Enrx(n,r,x): the weight enumerator of 132-avoiding 1-1- mappings
#in {1..n}->{{1+r,..,n+r} according to number of fixed pts
#according to the simplified recurrence
Enrx:=proc(n,r,x)
local gu,k:
option remember:
if n=0 then
 RETURN(1):
fi:
if n<0 then
 RETURN(0):
fi:

gu:=0:

for k from 1 to r do
gu:=gu+Cn(k-1)*Enrx(n-k,r-k,x):
od:

for k from r+1 to trunc((n+r+1)/2) do
gu:=gu+Cn(k-1)*Enrx(n-k,k-r,x):
od:

for k from trunc((n+r+1)/2)+1 to n do
gu:=gu+Cn(n-k)*Enrx(k-1,n-k+r,x):
od:
gu:=expand(gu):

if r=0 then
gu:=expand(gu+(x-1)*Enrx(n-1,0,x)):
fi:

gu:
end:

enrx:=proc(n,r,x)
local gu,i:
gu:=anx(n,x):
for i from 1 to r do
gu:=gu+(1-x)*Cn(i-1)*anx(n-i,x):
od:
expand(gu):
end:

#n is the length of the permutation
#k is number of fixed points
#w is number of occurrences
#pat from {abc,acb,bac,bca,cab,cba<->123,132,213,231,312,321}
numbrrp:=proc(pat,n,k,w)
local i,j,derset,dercount,S,T,flg,t,fixcount,D:
with(combinat):
S:=permute(n):T:={}:dercount:=0:D:={}:
for i from 1 to n! do
t:=S[i]: flg:=0:fixcount:=0:
for j from 1 to n do 
if t[j]=j then fixcount:=fixcount+1: fi:
if fixcount>k then flg:=1:fi:
od:
if fixcount<k then flg:=1:fi:
if flg=0 then T:=T union {t}:fi:
od:
for i from 1 to nops(T) do
t:=T[i]:
if pat(t,n,w)=1 then
dercount:=dercount+1:D:=D union {t}
fi:
od:
#D:
dercount:
end:

#123 Pattern
abc:=proc(prm,n,w)
local i,T,t,cnt:
T:=choose([seq(i,i=1..n)],3): cnt:=0:
for i from 1 to nops(T) do
t:=T[i]:
if (prm[t[1]]<prm[t[2]] and prm[t[2]]<prm[t[3]]) then
cnt:=cnt+1:
fi:
if cnt>w then RETURN(0):fi:
od:
if cnt=w then 1: else 0:fi:
end:

#132 Pattern
acb:=proc(prm,n,w)
local i,T,t,cnt:
T:=choose([seq(i,i=1..n)],3):cnt:=0:
for i from 1 to nops(T) do
t:=T[i]:
if (prm[t[1]]<prm[t[3]] and prm[t[3]]<prm[t[2]]) then 
cnt:=cnt+1:
if cnt>w then RETURN(0):fi:
fi:od:
if cnt=w then 1: else 0:fi:
end:
 
#213 Pattern
bac:=proc(prm,n,w)
local i,T,t,cnt:
T:=choose([seq(i,i=1..n)],3):cnt:=0:
for i from 1 to nops(T) do
t:=T[i]:
if (prm[t[2]]<prm[t[1]] and prm[t[1]]<prm[t[3]]) then
cnt:=cnt+1:
if cnt>w then RETURN(0):fi:
fi:od:
if cnt=w then 1: else 0:fi:
end:

#231 Pattern
bca:=proc(prm,n,w)
local i,T,t,cnt:
T:=choose([seq(i,i=1..n)],3):cnt:=0:
for i from 1 to nops(T) do
t:=T[i]:
if (prm[t[3]]<prm[t[1]] and prm[t[1]]<prm[t[2]]) then
cnt:=cnt+1:
if cnt>w then RETURN(0):fi:
fi:od:
if cnt=w then 1: else 0:fi:
end:
 
#312 Pattern
cab:=proc(prm,n,w)
local i,T,t,cnt:
T:=choose([seq(i,i=1..n)],3):cnt:=0:
for i from 1 to nops(T) do
t:=T[i]:
if (prm[t[2]]<prm[t[3]] and prm[t[3]]<prm[t[1]]) then
cnt:=cnt+1:
if cnt>w then RETURN(0):fi:
fi:od:
if cnt=w then 1: else 0:fi:
end:
 
#321 Pattern
cba:=proc(prm,n,w)
local i,T,t,cnt:
T:=choose([seq(i,i=1..n)],3):cnt:=0:
for i from 1 to nops(T) do
t:=T[i]:
if (prm[t[3]]<prm[t[2]] and prm[t[2]]<prm[t[1]]) then
cnt:=cnt+1:
if cnt>w then RETURN(0):fi:
fi:od:
if cnt=w then 1: else 0:fi:
end:

#Values of b(n,k) as defined in Lemma 4.2
bnk:=proc(n,k)
local i,j: 
with(combinat):
if k>n then RETURN(0):fi:
if k=n then RETURN(1):fi: 
if k=0 then RETURN(Cn(n)):fi:
bnk(n,k+1)+bnk(n-1,k-1):
end:

#Check for Lemma 4.2
bnkchk:=proc(n) 
local k:
add(bnk(n,k)*(x-1)^k,k=0..n):
sort(expand(%,x)):
end:


#Check for Lemma 4.3
bnk2:=proc(n,k)
local i:
with(combinat):
if k=0 then RETURN(Cn(n)):fi:
add(bnk(n-i,k-1)*Cn(i-1))/i,i=1..n):
end:

#Check for Lemma 7.4
bnk3:=proc(n,k)
with(combinat):
(k+1)/(n+1)*numbcomb(2*n-k,n):
end:

#Check for Cor 7.6
catchk:=proc(n)
local k,j:
add(add((-1)^j*numbcomb(j+k,k)*bnk3(n,k+j),j=0..n-k),k=0..n):
end:

#Check for Cor 7.7
fnchk:=proc(n)
local j:
add((-1)^j*bnk3(n,j),j=0..n):
end: