#HCRV.txt: One of the Maple programs accompanying the article
#"Counting Permutations that Avoid Many Patterns" by
#Yonah BIERS-ARIEL, Haripriya CHAKRABORTY, John CHIARELLI, Bryan EK, Andrew LOHR, 
#Jinyoung PARK, Justin SEMONSEN, Richard VOEPEL, Mingjia YANG, Anthony ZALESKI, 
#and Doron ZEILBERGER 
#Available from:
#http://www.math.rutgers.edu/~zeilberg/mamarim/mamarimhtml/pamp.html

with(combinat):
with(numtheory):

HelpHV:=proc():
print(``):
end:


#redu(L): For a list of reals L, outputs the reduction as a permutation of nops(L).
redu:=proc(L) local L1,T,k,i:
L1:=sort(sort(L,'output'='permutation'),'output'='permutation'):
end:


#IV(n,k): The set of increasing vectors of length k in {1, ...,n} ending in n.
IV:=proc(n,k) local S,i,s:
S:=choose([seq(i,i=1..n-1)],k-1):
{seq([op(s), n], s in S)}:
end:


#IsBad1(pi,p): Inputs a perm. pi and a pattern p, outputs true iff pi contains the
#pattern p where the last entry of pi participates.
IsBad1:=proc(pi,p) local S,n,k,s:
n:=nops(pi):
k:=nops(p):
S:=IV(n,k):
for s in S do
if redu(pi[s])=p then
RETURN(true):
fi:
od:
false: 
end:


#Inputs a non-neg. integer n and a set of patterns P and outputs the
#set of permutations avoiding the patterns of P 
GG:=proc(n,P) local Sold,S, i,cand,pi:
option remember:
if n=0 then
RETURN({[]}):
fi:
Sold:=GG(n-1,P):
S:={}:
for pi in Sold do
for i from 0 to n-1 do
cand:=redu([op(pi),i+1/2]):
if not IsBad(cand,P) then
S:=S union {cand}:
fi:
od:
od:
S:
end:


#SeqGG(P,N): Inputs a set of patterns P, and a positive integer N
#returns the list of length N whose i-th entry is the NUMBER of permutations
#of length i NOT containing any of the patterns in P
SeqGG:=proc(P,N) local i:
[seq( nops(GG(i,P)),i=1..N)]:
end:


#GP1(L,n,d): Inputs a list of pairs [input,output], and a variable n
#and outputs a polynomial of degree d such P(input)=output
GP1:=proc(L,n,d) local P,i,a,var,eq:
if d>=nops(L)-2 then
RETURN(FAIL):
fi:
P:=add(a[i]*n^i,i=0..d):
var:={seq(a[i],i=0..d)}:
eq:={seq(subs(n=L[i][1],P)=L[i][2],i=1..d+2)}:
var:=solve(eq,var):
if var=NULL then
RETURN(FAIL):
fi:
P:=subs(var,P):
if {seq(subs(n=L[i][1],P)-L[i][2],i=d+3..nops(L))}<>{0} then
RETURN(FAIL):
else
RETURN(P):
fi:
end:


#GP(L,n): inputs a list of pairs [input,output], and a variable n
#and outputs a polynomial of degree <=nops(L) such P(input)=output
GP:=proc(L,n) local d,hope:
for d from 0 to nops(L)-2 do
hope:=GP1(L,n,d):
if hope<>FAIL then
RETURN(hope):
fi:
od:
FAIL:
end:


#Inputs a permutation pi and a set of patterns P and outputs true if
#it contains at least one member of P.
IsBad:=proc(pi,P) local p:
for p in P do
if IsBad1(pi,p) then
RETURN(true):
fi:
od:
false:
end:


###########################################################################


#revers(perm): The reverse of a permutation.
revers:=proc(perm) local i:
[seq(perm[nops(perm)-i+1],i=1..nops(perm))]:
end:


#REVERS(Perms): The set of reverses of a set of permutations Perms.
REVERS:=proc(Perms) local i,gu:
gu:={}:
for i from 1 to nops(Perms) do
gu:=gu  union{revers(op(i,Perms))}:
od:
gu:
end:


hofkhi:=proc(pi) local i,T:
for  i from 1 to nops(pi) do
T[pi[i]]:=i:
od:
[seq(T[i],i=1..nops(pi))]:
end:


Hofkhi:=proc(S) local i:
{seq(hofkhi(S[i]),i=1..nops(S))}:
end:


INVERSE:=proc(Perms) local i:
{seq(hofkhi(Perms[i]),i=1..nops(Perms))}:
end:


#KHAVERIM(Perms): Given a set of patterns Perms, outputs all the images under the 
#trivial Wilf-equivalence group. For example, try: KHAVERIM({[1,2,3]});
KHAVERIM:=proc(Perms):
{Perms, REVERS(Perms),INVERSE(Perms),
REVERS(INVERSE(Perms)),INVERSE(REVERS(Perms)),
INVERSE(REVERS(INVERSE(Perms))),
REVERS(INVERSE(REVERS(Perms))),
INVERSE(REVERS(INVERSE(REVERS(Perms))))}:
end:


#MakeRPS(m,n): Constructs a set containing representatives from each  
#trivial Wilf-equivalence class on sets of m permutations on n elements.
MakeRPS:=proc(m,n) local S,RPS,bool,s,k:
option remember:
S:=choose({op(permute(n))},m):
RPS:={}:
bool:=false:
for s in S do
for k in KHAVERIM(s) do
if k in RPS then
bool:=true:
fi:
od:
if not bool then
RPS:= RPS union {s}:
fi:
bool:=false:
od:
RETURN(RPS):
end:


#Chunk(S): Takes a set S and splits it into a number of dijoint subsets of size
#at most floor(sqrt(nops(S))), returning an indexed collection SC[i]. For example,
#a set of size 14 is chunked into 4 sets of size 3 and one set of size 2, 
#and the output SC has 4 defined indices.
Chunk:=proc(S) local n,T,SC,i,t:
n:=floor(sqrt(nops(S))):
T:=S:
i:=1:
while nops(T)>=n do
SC[i]:={}:
while nops(SC[i])<n do
t:=T[1]:
SC[i]:=SC[i] union {t}:
T:=T minus {t}:
od:
i:=i+1:
od:
SC[i]:=T:
RETURN(SC):
end: