######################################################################
##P12345: Save this file as  P12345                                  #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read P12345<Enter>                                   #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Rutgers University ,                  #
#zeilberg at math dot rutgers dot edu                                #
######################################################################
 
#Created: June-Aug. 2012
with(combinat):
 
print(`Created:  June-Aug. 2012`):
print(` This is P12345, one of the numerous Maple packages `):
print(`accompanying Brian Nakamura and Doron Zeilberger's  article: `):
print(`"Using Noonan-Zeilberger Functional Equations `):
print(`to enumerate (in Polynomial Time!) Generalized Wilf classes`):
print(`available from:`):
print(` http://www.math.rutgers.edu/~zeilberg/mamarim/mamarimhtml/Gwilf.html .`):
print():
print(`Please report bugs to zeilberg at math dot rutgers dot edu`):
print(``):
 print(`The most current version of this  package `):
 print(` is  available from`):
 print(`http://www.math.rutgers.edu/~zeilberg/tokhniot/P12345  .`):



print(`-------------------------------------------------------------`):
 print(`For a list of the Checking procedures type ezraC();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`-------------------------------------------------------------`):

print(`-------------------------------------------------------------`):
 print(`For a list of the supporting procedures type ezra1();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`-------------------------------------------------------------`):


print(`-------------------------------------------------------------`):
 print(`For a list of the MAIN procedures type ezra();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`-------------------------------------------------------------`):

ezraC:=proc()

if args=NULL then
 print(` The checking procedures are: CheckPn, PnDirect `):

else
ezra(args):
fi:

end:




ezra1:=proc()

if args=NULL then
 print(` The supporting procedures are: Children `):
 print(`  LtoS  `):
 print(`  PqLyT, TrimTriple, TrimTriples, Wt, Yeladim, YeladimT `):

else
ezra(args):
fi:

end:

ezra:=proc()

if args=NULL then
 print(`The main procedures are:`):
 print(` fn, `):
 print(` Pn, Seqs `):
 print(` `):



elif nops([args])=1 and op(1,[args])=CheckPn then
print(`CheckPn(N):Checks Pn for n from 1 to N `):
print(`Try: `):
print(` CheckPn(8); `):



elif nops([args])=1 and op(1,[args])=fn then
print(`fn(n,q): the weight-enumerator of S_n according to the weight`):
print(`q^(#12345)`):


elif nops([args])=1 and op(1,[args])=LtoS then
print(`LtoS(L,q,r): given a vector of increasing powers of q up to the`):
print(`q^(r+1) returns the list of integers of length r+2 whose`):
print(`j-th entry is the number of times q^j shows up`):
print(`Try LtoS([1,1,q,q,q^3,q^3,q^3],q,2);`):




elif nops([args])=1 and op(1,[args])=Pn then
print(`Pn(n,q,z1,z2,z3):The weight-enumerator of perms of {1, ..., n}`):
print(`according to the weight Wt(pi,q,z1,z2,z3) ...`):
print(`try: `):
print(`Pn(5,q,z1,z2,z3); `):

elif nops([args])=1 and op(1,[args])=PnDirect then
print(`PnDirect(n,q,z1,z2,z3):The weight-enumerator of perms of {1, ..., n}`):
print(`according to the weight Wt(n,q,z1,z2,z3) `):
print(`computed directly. FOR CHECKING PURPOSES ONLY, VERY SLOW!`):
print(`try: `):
print(`PnDirect(5,q,z1,z2,z3); `):

elif nops([args])=1 and op(1,[args])=PqL then
print(` PqL(q,L1,L2):The value of P_n(n,q,y,z) where n=nops(L1)=nops(L2)`):
print(`and one plugs in for [y[1], ..., y[n]] the`):
print(`values of L1, and for [z[1], ..., z[n]] the values of L2 try:`):
print(`PqL(q,[y1,y2,y3,y4],[z1,z2,z3,z4]);`):


elif nops([args])=1 and op(1,[args])=PqLyT then
print(`PqLyT(q,L1,L2,L3,r): only interested in powers up to r`):
print(`to the r-th power of q`):
print(`Try:`):
print(`PqLyT(q,[1,1,1,1],[1,1,1,1],[1,1,1,1],1);`):


elif nops([args])=1 and op(1,[args])=Seqs then
print(`Seqs(N,R): The lists whose R+1 lists whose (r+1)-th list`):
print(`is the first N members (starting at n=1)`):
print(`number of permutations with exactly r occurrence of 12345,`):
print(`try: Seqs(8,3);`):




elif nops([args])=1 and op(1,[args])=TrimTriple then
print(`TrimTrible(P,q,r): Given a Triple P=[coeff,L1,L2]`):
print(`where coeff is a power of q, `):
print(` and L1, L2 are lists of powers of q, of the same length`):
print(`returns the `):
print(`empty set if coeff has power larger than r,`):
print(`and otherwise returns the singleton`):
print(`{[coeff,L1a,L2a]} where L1a is L1 with`):
print(`any power larger than q^(r+1) is replaced by q^(r+1)`):
print(`and ditto for L2a and L2`):
print(`For example, try: `):
print(`TrimTriple([q,[1,1,q,q,q^3,q^4,q^4],[1$7]],q,2);`):

elif nops([args])=1 and op(1,[args])=TrimTriples then
print(`TrimTriples(S,q,r): the set of trimmed triples coming from the`):
print(`set of triples S with respect to the power r, try:`):
print(`TrimTriples({[q,[1,1,q,q,q^3,q^4,q^4]]},q,1);`):

elif nops([args])=1 and op(1,[args])=Wt then
print(`Wt(pi,q,z1,z2,z3): the weight of a permutation pi according to`):
print(`q^(#12345)*`):
print(`Product(z1[i]^#(1234,1=i)*z2[i]^#(123,1=i) ... `):
print(`*z3[i]^#(12,i=1), i=1..n) `):
print(`Try: `):
print(`Wt([1,2,3,4,5],q,z1,z2,z3); `):



elif nops([args])=1 and op(1,[args])=Yeladim then
print(`Yeladim(L1,L2,L3,q): the children of the lists L1,L2,L3`):
print(`Try: Yeladim([1,1,1],[1,1,1],[1,1,1],q);`):

elif nops([args])=1 and op(1,[args])=YeladimT then
print(`YeladimT(L1,L2,L3,q,r): the trimmed children of the list L1,L2,L3`):
print(`with respect to the power r.`):
print(`Try: YeladimT([1,q^3,q^5],[1,q,q^3],[1,1,1],q,2);`):

else
print(`There is no ezra for`,args):
fi:
 
end:


#fn(n,q): the weight-enumerator of S_n according to the weight
#q^(#12345)
fn:=proc(n,q) local gu,i,z1,z2,z3:
gu:=Pn(n,q,z1,z2,z3):
sort(subs({seq(z1[i]=1,i=1..n), seq(z2[i]=1,i=1..n)
, seq(z3[i]=1,i=1..n)},gu)):
end:

#Pn(n,q,z1,z2,z3):The weight-enumerator of perms of {1, ..., n}
#according to the weight Wt(pi,q,z1,z2,z3)
Pn:=proc(n,q,z1,z2,z3) local i,j,gu:
option remember:
if n=1 then
 RETURN(1):
else
 gu:=Pn(n-1,q,z1,z2,z3):
 expand(add(z3[i]^(n-i)*subs(
{
seq(z1[j]=q*z1[j+1],j=i..n-1),
seq(z2[j]=z1[i]*z2[j+1],j=i..n-1),
seq(z3[j]=z2[i]*z3[j+1],j=i..n-1)
},gu),i=1..n)):
fi:
end:


#Wt(pi,q,z1,z2,z3): the weight of a permutation pi according to
#q^(#12345)*z1[i]^#(1234,1=i)*z[2][i]^#(123,1=i) ... `):
#*z3[i][1]#(12,1=i)
#Try: 
#Wt([1,2,3,4,5],q,z1,z2,z3);
Wt:=proc(pi,q,z1,z2,z3) local n,i1,i2,i3,i4,i5,mu:
n:=nops(pi):

mu:=1:

for i1 from 1 to n do
 for i2 from i1+1 to n do
   for i3 from i2+1 to n do
    for i4 from i3+1 to n do
     for i5 from i4+1 to n do
    if pi[i1]<pi[i2] and pi[i2]<pi[i3] and pi[i3]<pi[i4] and pi[i4]<pi[i5] 
       then
          mu:=mu*q:
    fi:
   od:
 od:
od:
od:
od:

for i1 from 1 to n do
 for i2 from i1+1 to n do
  for i3 from i2+1 to n do
   for i4 from i3+1 to n do
   if pi[i1]<pi[i2] and pi[i2]<pi[i3] and pi[i3]<pi[i4] then
     mu:=mu*z1[pi[i1]]:
   fi:
 od:
od:
od:
od:


for i1 from 1 to n do
 for i2 from i1+1 to n do
  for i3 from i2+1 to n do
   if pi[i1]<pi[i2] and pi[i2]<pi[i3] then
     mu:=mu*z2[pi[i1]]:
   fi:
 od:
od:
od:

for i1 from 1 to n do
 for i2 from i1+1 to n do
   if pi[i1]<pi[i2] then
     mu:=mu*z3[pi[i1]]:
   fi:
od:
od:

mu:
end:


#PnDirect(n,q,z1,z2,z3):The weight-enumerator of all perms of
#of length n
#according to Wt
#done directly, for checking purposes only! Very slow
PnDirect:=proc(n,q,z1,z2,z3) local gu,pi:
gu:=permute(n):
add(Wt(pi,q,z1,z2,z3), pi in gu):
end:



#CheckPn(N):Checks Pn(n,q,y) for n from 1 to N 
#Try:
#CheckPn(8);
CheckPn:=proc(N) local n,q,z1,z2,z3:
 evalb({seq(PnDirect(n,q,z1,z2,z3)-Pn(n,q,z1,z2,z3),n=1..N)}={0}):
end:





#Seqs(N,R): The lists whose R+1 lists whose (r+1)-th list
#is the first N members (starting at n=1)
#number of permutations with exactly r occurrence of 12345,
#try: Seqs(8,3);
Seqs:=proc(N,R) local gu,n1,q,r:
option remember:
gu:=[seq(PqLyT(q,[1$n1],[1$n1],[1$n1],R),n1=1..N)]:
[seq([seq(coeff(gu[n1],q,r),n1=1..N)],r=0..R)]:
end:

#PqLyT(q,L1,L2,L3,r):Like PqLy but only interested in terms up
#to the r-th power of q
#Try:
#PqLyT(q,[1,1,1,1],[1,1,1,1],[1,1,1,1,1],1);
PqLyT:=proc(q,L1,L2,L3,r) local n,i,gu,mu:
option remember:
n:=nops(L1):
if nops(L2)<>n or nops(L3)<>n then
 print(`Bad input`):
 RETURN(FAIL):
fi:

if n=1 then
 RETURN(1):
else

gu:=YeladimT(L1,L2,L3,q,r):
mu:=expand(add(gu[i][1]*PqLyT(q,gu[i][2],gu[i][3],gu[i][4],r),i=1..nops(gu))):
mu:=add(coeff(mu,q,i)*q^i,i=0..r):
RETURN(mu):
fi:
end:

#TrimTriple(P,q,r): Given a Quad P=[coeff,L1,L2,L3]
#where coeff is a power of q, 
#and L1, L2, L3 are lists of powers of q, returns the
#empty set if coeff has power larger than r,
#and otherwise returns the singleton
#{[coeff,L1a,L2a,L3a]} where L1a is L1 with
#any power larger than q^(r+1) is replaced by q^(r+1)
#and ditto for L2, and L3
#For example, try: 
#TrimTriple([q,[1,1,q,q,q^3,q^4,q^4],[1,1,q,q^2,q^4,q^5],[1$7],q,2);
TrimTriple:=proc(P,q,r) local c,L1,L2,L3,L1a,L2a,L3a,i:
c:=P[1]:
L1:=P[2]:
L2:=P[3]:
L3:=P[4]:
if degree(c,q)>r then
 RETURN({}):
fi:

L1a:=[]:

for i from 1 to nops(L1) do
 if degree(L1[i],q)>r+1 then
  L1a:=[op(L1a),q^(r+1)]:
 else
  L1a:=[op(L1a),L1[i]]:
 fi:
od:


L2a:=[]:
for i from 1 to nops(L2) do
 if degree(L2[i],q)>r+1 then
  L2a:=[op(L2a),q^(r+1)]:
 else
  L2a:=[op(L2a),L2[i]]:
 fi:
od:

L3a:=[]:
for i from 1 to nops(L3) do
 if degree(L3[i],q)>r+1 then
  L3a:=[op(L3a),q^(r+1)]:
 else
  L3a:=[op(L3a),L3[i]]:
 fi:
od:


{[c,L1a,L2a,L3a]}:
end:



#TrimTriples(S,q,r): the set of trimmed pairs coming from the
#set of pairs S with respect to the power r, try:
#TrimTripless({[q,[1,1,q,q,q^3,q^4,q^4]]},q,1);
TrimTriples:=proc(S,q,r) local P:
{seq(op(TrimTriple(P,q,r)), P in S)}:
end:

#Yeladim(L1,L2,L3,q): the children of the lists L1,L2,L3
#Try: Yeladim([1,1,1],[1,1,1],[1,1,1],q);
Yeladim:=proc(L1,L2,L3,q) local n,i,j:
n:=nops(L1):
if n<>nops(L2) then
 print(`Bad input`):
 RETURN(FAIL):
fi:

[seq([L3[i]^(n-i),[op(1..i-1,L1),seq(q*L1[j],j=i+1..n)],
[op(1..i-1,L2),seq(L1[i]*L2[j],j=i+1..n)],
[op(1..i-1,L3),seq(L2[i]*L3[j],j=i+1..n)]
],i=1..n)];
end:


#YeladimT(L1,L2,L3,q,r): the trimmed children of the lists L1,L2,L3
#with respect to the power r.
#Try: YeladimT([1,q^3,q^5],[1,q,q^2],[1,q,q^2], q,2);
YeladimT:=proc(L1,L2,L3,q,r):
TrimTriples(Yeladim(L1,L2,L3,q),q,r):
end: