######################################################################
## NewWilf.txt Save this file as  NewWilf.txt  to use it,            #
# stay in the                                                        #
## same directory, get into Maple (by typing: maple <Enter> )        #
## and then type:  read  `NewWilf.txt` <Enter>                       #
## Then follow the instructions given there                          #
##                                                                   #
## Written by Lucy Martinez and Doron Zeilberger, Rutgers University #
######################################################################

with(combinat):

print(`First Written: May, 2025: tested for Maple 2020 `):
print(`Version: May 22, 2025 `):
print():
print(`This is NewWilf.txt, A Maple package`):
print(`accompanying Lucy Martinez and Doron Zeilberger's  article: `):
print(` Counting Parking Functions Whose Permutations Avoid Patterns`):
print(`This package has nothing to do with Parking functions but it would`):
print(`enable to be generalized to parking functions `):
print(`It is one of the packages available from `):
print(` http://sites.math.rutgers.edu/~zeilberg/tokhniot/mamarim/mamarimhtml/pwilf.html `):
print():
print(`The most current version is available on WWW at:`):
print(` http://sites.math.rutgers.edu/~zeilberg/tokhniot/NewWilf.txt .`):
print(`Please report all bugs to: DoronZeil at gmail dot com .`):
print():

print(`-------------------------------------------`):
print(`For general help, and a list of the MAIN functions,`):
print(` type "ezra();". For specific help type "ezra(procedure_name);" `):
print(`For a list of the supporting functions type: ezra1();`):
print():
print(`-------------------`):


ezra1:=proc()
if args=NULL then

print(`The SUPPORTING procedures are: DomainWilf, FindGaps, FindGaps1,  Insert, IsGood, IsRevDel, IV, redu, RevDelSet, WilfSeqBF, WilfSeqVS, WilfTail `):
print(``):

else
ezra(args):

fi:

end:


ezra:=proc()
if args=NULL then

print(`NewWilf.txt: A Maple package for  counting, using Experimental Mathematics, Wilf classes of permuations.`):

print(`The MAIN procedures are: EvalScheme, FindScheme, SchemeSeq, Wilf`):



elif nargs=1 and args[1]=DomainWilf then
print(`DomainWilf(n,Patterns,sig): inputs a pos. integer n, a set of Patterns, Patterns, and a permutation sig, outputs the set of tails`):
print(`in Wilf(n,Patterns) that reduce to sig. Try:`):
print(`DomainWilf(7,{[1,2,3]},[1,2]);`):




elif nargs=1 and args[1]=FindGaps then
print(`FindGaps(n,Patterns,sig,r): Looks at FindGaps1(n1,Patterns,sig) for n1=n,n+1,n+2, and outputs a vector of length nops(sig)+1 where 1/2 means no agreement otherwise that they agree `):
print(`FindGaps(6,{[1,2,3],[2,3,1]},[2,1],2);`):



elif nargs=1 and args[1]=FindGaps1 then
print(`FindGaps1(n,Patterns,sig): Inputs a set of patterns, Patterns, a positive integer n, Looks at DomainWilf(n, Patterns,sig), and looks the`):
print(`set of (let k:=nops(sig)) length k+1 {seq(w[1]-1,[seq(w[i]-w[i-1],i=2..nops(w)),n-w[k]], w) for all w in DomainWilf(n,Patterns,sig)`):
print(`and takes the minumum of eacu entry. The output is a list of lenth k+1. Try:`):
print(`FindGaps1(8,{[1,2,3],[2,3,1]},[2,1]);`):



elif nargs=1 and args[1]=Insert then
print(`Insert(pi,i): Given a permutation pi, and an integer i between 1 and nops(pi)+1, constructs the permutation`):
print(`of length nops(pi)+1 whose last entry is i, and such that reduced form of the first nops(pi) letters is pi. Try`):
print(`Insert([1,4,2,3],2);`):
 
elif nargs=1 and args[1]=IsGood then
print(`IsGood(pi,SetPatterns): Given a permutation pi, and a set of patterns SetPatterns, decides whether`):
#pattern1 occurs in pi with the first and last letter of
print(`pi  coinciding with the last letter of pattern1. It returns false if this is the case, otherwise true`):

elif nargs=1 and args[1]=IsRevDel then
print(`IsRevDel(n,Patterns,sig,k): Is the k-th entry of all possible tails of length nops(sig) that reduce to sig Reversely deletable. Try:`):
print(`IsRevDel(8,{[1,2,3]},[2,1],2);`):

elif nargs=1 and args[1]=IV then
print(`IV(n,k) all increasing vectors of length k of the form 1<=i_1< ...<i_k <=n. Try:`):
print(`IV(4,2);`):


elif nargs=1 and args[1]=FindScheme then
print(`FindScheme(n,Patterns,k,r): Tries to find a scheme of length k. Tested for n,n+1,..,n+r. Try:`):
print(`FindScheme(7,{[1,2,3]},2,2); The output is a list of length 8`):
print(`[k,S1,op(T1),op(T3),S2,op(T4),n1, WilfSeqBF(Patterns,n1-1)]:`):
print(`where `):
print(`k is the size of the suffix, `):
print(`S1 is the set of reductions of the k-suffices of members of Wilf(n,Patterns)`):
print(`op(T1) is the reduction table`):
print(`op(T3) is the Gap vectors Table`):
print(`S2 is the set of suffixes for the initial case n=n1`):
print(`op(T4) is the accompanying table`):
print(`n1 is where the scheme starts to be valid`):
print(` WilfSeqBF(Patterns,n1-1): the first n1-1 terms of the enumerating sequence. Try:`):
print(`FindScheme(7,{[1,2,3]},2,1);`):


elif nargs=1 and args[1]=EvalScheme then
print(`EvalScheme(S,n,L): If a  scheme S=[k,S1,op(T1),op(T2),op(T3),S2,op(T4),n1]  is supposed to count the number of permutations avoiding the set of patterns Patterns`):
print(`then the output is the number of such permutations of length n whose last entries is L. Try:`):
print(`S:=FindScheme(7,{[1,2,3]},2,2); EvalScheme(S,8,[6,3]);`):


elif nargs=1 and args[1]=EvalSchemeA then
print(`EvalSchemeA(S,n): If a  scheme S=[k,S1,op(T1),op(T2),op(T3),S2,op(T4),n1]  is supposed to count the number of permutations avoiding the set of patterns Patterns`):
print(`then the output is the number of such permutations of length n. Try:`):
print(`S:=FindScheme(7,{[1,2,3]},2,2); EvalSchemeA(S,8);`):

elif nargs=1 and args[1]=redu then
print(`redu(perm): Given a permutation on a set of positive integers, finds its reduced form. Try:`):
print(`redu([5,7,3]);`):

elif nargs=1 and args[1]=RevDelSet then
print(`RevDelSet(n,Patterns,sig): all the k such that k is reversely deletable for sig. Try`):
print(`RevDelSet(6,{[1,2,3]},[2,1]);`):


elif nargs=1 and args[1]=SchemeSeq then
print(`SchemeSeq(S,N): the first N terms evaluated by the scheme S. Try:`):
print(`S:=FindScheme(7,{[1,2,3]},2,1); SchemeSeq(S,20);`):

elif nargs=1 and args[1]=Wilf then
print(`Wilf(n,Patterns): inputs a non-neg.  integer n and a set of Patterns, outputs the set if permutations of {1, ...,n} that avoid all the patterns in Patterns. Try: `):
print(`Wilf(7,{[1,2,3],[1,3,2]});`):

elif nargs=1 and args[1]=WilfTail then
print(`WilfTail(n,Patterns,zanav): all the members of Wilf(n,Patterns) that end in zanav. Try:`):
print(`WilfTail(7,{[3,1,2]},[1,2]);`):

elif nargs=1 and args[1]=WilfSeqBF then
print(`WilfSeqBF(Patterns,N): The first N terms, by brute force. Try:`):
print(`WilfSeqBF({[1,2,3]}, 8);`):


elif nargs=1 and args[1]=WilfSeqVS then
print(`WilfSeqVS(n,Patterns,k,r,N): The first N terms, by brute force. Try:`):
print(`WilfSeqVS(7,{[1,2,3]},2,1,8): First tries to find a scheme via FindScheme(n,Patterns,k,r) and then uses it to compute the first N terms`):

else
 print(`There is no such thing as`, args):

fi:

end:


##start from Wilf.txt 

#IV(n,k) all increasing vectors of length
#k of the form 1<=i_1< ...<i_k <=n

IV:=proc(n,k)
local i,gu,mu:

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

if n<k then
RETURN({}):
fi:

gu:=IV(n-1,k):
mu:=IV(n-1,k-1):

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

gu:
end:


#IsGood1(pi,pattern1)
#Given a permutation pi, and a pattern pattern1, decides whether
#pattern1 occurs in pi with the first and last letter of pi  coinciding
#with the last letter of pattern1. It returns false if this is
#the case, otherwise true

IsGood1:=proc(pi,pattern1)
local gu,k,n,i,j,subperm,vec:

if pi=pattern1 then
RETURN(0):
fi:

n:=nops(pi):

k:=nops(pattern1):
if k>n then
RETURN(true):
fi:
gu:=IV(n,k-2):

for i from 1 to nops(gu) do
vec:=op(i,gu):
subperm:=[op(1,pi)]:
for j from 1 to k-2 do
subperm:=[op(subperm),pi[vec[j]]]:
od:
subperm:=[op(subperm),pi[n]]:
if redu(subperm)=pattern1 then
RETURN(false):
fi:
od:

true:

end:



#IsGood(pi,SetPatterns)
#Given a permutation pi, and a set of patterns
#SetPatterns, decides whether
#pattern1 occurs in pi with the first and last letter of
#pi  coinciding
#with the last letter of pattern1. It returns false if this is
#the case, otherwise true

IsGood:=proc(pi,SetPatterns) local i:
if member(pi,SetPatterns) then
RETURN(false):
fi:

for i from 1 to nops(SetPatterns) do
if not IsGood1(pi,op(i,SetPatterns)) then
 RETURN(false):
fi:
od:
true:
end:
 
 
#redu(perm): Given a permutation on a set of positive integers
#finds its reduced form
 
redu:=proc(perm)
local gu,i,ra,mu:
gu:=sort(perm):
 
for i from 1 to nops(gu) do
ra[gu[i]]:=i:
od:
 
mu:=[]:
 
for i from 1 to nops(perm) do
mu:=[op(mu),ra[perm[i]]]:
od:
mu:
end:
 
#Insert(pi,i): Given a permutation pi, and an integer
#i between 1 and nops(pi)+1, constructs the permutation
#of length nops(pi)+1 whose last entry is i, and such that
#reduced form of the first nops(pi) letters is pi
 
Insert:=proc(pi,i)
local mu,j:
mu:=[]:
 
for j from 1 to nops(pi) do
if pi[j]<i then
mu:=[op(mu),pi[j]]:
else
mu:=[op(mu),pi[j]+1]:
fi:
od:
mu:=[op(mu),i]:
mu:
end:

#Wilf(n,Patterns): all permutations of length n
#that avoid the patterns in the set of patterns
#Patterns
 
Wilf:=proc(n,Patterns)
local i,mu,gu,pi,j,muamad:
option remember:
 
 
if n=0 then
RETURN({[]}):
fi:
 
mu:=Wilf(n-1,Patterns):
gu:={}:
 
for i from 1 to nops(mu) do
pi:=op(i,mu):
 
for j from 1 to n do
muamad:=Insert(pi,j):
if member(redu([op(2..n,muamad)]),mu) then
  if IsGood(muamad,Patterns) then
   gu:=gu union {muamad}:
  fi:
fi:
od:
od:
 
gu:
end:




##End from Wilf.txt 


#WilfTailTable(n,Patterns,r): all the tails of size r of Wilf(n,Patterns), followed by the table. Try:
#WilfTailTable(7,{[3,1,2]},2);
WilfTailTable:=proc(n,Patterns,r) local gu,gu1,mu,mu1,T:
option remember:
mu:=Wilf(n,Patterns):
if r>n then
 RETURN(FAIL):
fi:

gu:={seq([op(n-r+1..n,mu1)], mu1 in mu)}:

for gu1 in gu do
 T[gu1]:={}:
od:

for mu1 in mu do
T[[op(n-r+1..n,mu1)]]:=T[[op(n-r+1..n,mu1)]] union {mu1}:
od:
gu,op(T):

end:

#WilfTail(n,Patterns,zanav): all the members of Wilf(n,Patterns) that end in zanav. Try:
#WilfTail(7,{[3,1,2]},[1,2]);
WilfTail:=proc(n,Patterns,zanav) local r,gu:
r:=nops(zanav):
gu:=WilfTailTable(n,Patterns,r):
if not member(zanav,gu[1]) then
 {}:
else
gu[2][zanav]:
fi:
end:

#DomainWilfTable(n,Patterns,r): The table that for each permutation of r gives you all tails of length r of members of Wilf(n,Patterns) that reduce to it. Try:
#DomainWilfTable(8,{[1,2,3]},2);
DomainWilfTable:=proc(n,Patterns,r) local gu,gu1,T,mu,mu1:
option remember:

gu:=permute(r):

for gu1 in gu do
 T[gu1]:={}:
od:

mu:=WilfTailTable(n,Patterns,r)[1]:

for mu1 in mu do
 T[redu(mu1)]:= T[redu(mu1)] union {mu1}:
od:
op(T):
end:



#DomainWilf(n,Patterns,sig): inputs a pos. integer n, a set of Patterns, Patterns, and a permutation sig, outputs the set of tails
#in Wilf(n,Patterns) that reduce to sig. Try:
#DomainWilf(7,{[1,2,3]},[1,2]);
DomainWilf:=proc(n,Patterns,sig)
option remember:
DomainWilfTable(n,Patterns,nops(sig))[sig]:
end:


#IsRevDel1(n,Patterns,zanav,k): Is the k-th entry of zanav reversely deletable for permutations of length n that avoid the set of patterns Patterns and that end in zanav? Try
#IsRevDel1(8,{[2,1,3]},[5,2],2);
IsRevDel1:=proc(n,Patterns,zanav,k) local gu,z,gu1,lu,akha:
gu:=WilfTail(n,Patterns,zanav):
z:=nops(zanav):
if gu=FAIL then
  RETURN(FAIL):
fi:

if not (k>=1 and k<=z) then
 RETURN(FAIL):
fi:
gu:={seq(redu([op(1..n-z+k-1,gu1),op(n-z+k+1..n,gu1)]), gu1 in gu)}:

if z=1 then
lu:=Wilf(n-1,Patterns):
RETURN(evalb(gu=lu)):
else

akha:={seq([op(n-z+1..n-1,gu1)],gu1 in gu)}:
if nops(akha)<>1 then
 RETURN(FAIL):
fi:
akha:=akha[1]:
lu:=WilfTail(n-1,Patterns,akha):
evalb(gu=lu):
fi:
end:


#IsRevDel(n,Patterns,sig,k): Is the k-th entry of all possible tails of length nops(sig) that reduce to sig Reversely deletable. Try:
#IsRevDel(8,{[1,2,3]},[2,1],2);
IsRevDel:=proc(n,Patterns,sig,k) local gu,gu1:

gu:=DomainWilf(n, Patterns, sig):
evalb({seq(IsRevDel1(n,Patterns,gu1,k),gu1 in gu)}={true}):
end:



#FindGaps1(n,Patterns,sig): Inputs a set of patterns, Patterns, a positive integer n, Looks at DomainWilf(n, Patterns,sig), and looks the
#set of (let k:=nops(sig)) length k+1 {seq(w[1]-1,[seq(w[i]-w[i-1],i=2..nops(w)),n-w[k]], w) for all w in DomainWilf(n,Patterns,sig)
#and takes the minumum of eacu entry. The output is a list of lenth k+1. Try:
#FindGaps1(8,{[1,2,3],[2,3,1]},[2,1]);
FindGaps1:=proc(n,Patterns,sig) local gu,w,i,mu,k,mu1,gu1:
k:=nops(sig):

gu:=DomainWilf(n,Patterns,sig):

gu:={seq(sort(gu1),gu1 in gu)}:

mu:={seq([w[1]-1,seq(w[i]-w[i-1]-1,i=2..k), n-w[k]], w in gu)}:

[seq(max(seq(mu1[i],mu1 in mu)),i=1..k+1)]:

end:

#FindGaps(n,Patterns,sig,r): Looks at FindGaps1(n1,Patterns,sig) for n1=n,n+1,n+2, and outputs those that are the same for all. 
#FindGaps(6,{[1,2,3],[2,3,1]},[2,1],2);
FindGaps:=proc(n,Patterns,sig,r) local k,gu,n1,C,ku,gu1,i:

k:=nops(sig):

gu:={seq(FindGaps1(n1,Patterns,sig),n1=n..n+r)}:

C:=[];

for i from 1 to k+1 do
 ku:={seq(gu1[i],gu1 in gu)}:
 if nops(ku)=1 then
  C:=[op(C),ku[1]]:
 else
  C:=[op(C),1/2]:
 fi:
od:


end:




#RevDelSet1(n,Patterns,sig): all the k such that k is reversely deletable for sig
RevDelSet1:=proc(n,Patterns,sig) local k,gu:
gu:={}:
for k from 1 to nops(sig)  do
if IsRevDel(n,Patterns,sig,k) then
 gu:=gu union {k}:
fi:
od:
gu:
end:



#RevDelSet(n,Patterns,sig,r): all the k such that k is reversely deletable for sig for n to n+r
RevDelSet:=proc(n,Patterns,sig,r)  local n1,lu:
lu:={seq(RevDelSet1(n1,Patterns,sig),n1=n..n+r)}:
if nops(lu)<>1 then
 RETURN(FAIL):
fi:
for n1 from n by -1 while RevDelSet1(n1,Patterns,sig)=lu[1] do od:
lu[1],n1+1:
end:




#FindScheme(n,Patterns,k,r): Tries to find a scheme of length k. Tested for n,n+1,..,n+r. Try:
#FindScheme(7,{[1,2,3]},2,2); The output is a list of length 8
#[k,S1,op(T1),op(T3),S2,op(T4),n1, WilfSeqBF(Patterns,n1-1)]:
#where 
#k is the size of the suffix, 
#S1 is the set of reductions of the k-suffices of members of Wilf(n,Patterns)
#op(T1) is the reduction table
#op(T3) is the Gap vectors Table
#S2 is the set of suffixes for the initial case n=n1
#op(T4) is the accompanying table
#n1 is where the scheme starts to be valid
# WilfSeqBF(Patterns,n1-1): the first n1-1 terms of the enumerating sequence. Try:
#FindScheme(7,{[1,2,3]},2,1);

FindScheme:=proc(n,Patterns,k,r) local gu,S1,T1,T3,gu1,ka,gadol,sig,n1,S2,w,T4,S:
option remember:
gu:=Wilf(n+r,Patterns):
S1:={seq(redu([op(nops(gu1)-k+1..nops(gu1),gu1)]),gu1 in gu)}:
gu:=Wilf(n,Patterns):
if S1<>{seq(redu([op(nops(gu1)-k+1..nops(gu1),gu1)]),gu1 in gu)} then
 RETURN(FAIL):
fi:

gadol:=1:


for gu1 in S1 do
   ka:=RevDelSet(n,Patterns,gu1,r):
 if ka={} or ka[1]={} then
  RETURN(FAIL):
 else
 T1[gu1]:=ka[1]:
   gadol:=max(gadol,ka[2]):
 fi:
od:


for sig in S1 do
ka:=FindGaps(n,Patterns,sig,r):
T3[sig]:=ka:
od:




S2:={}:
n1:=max(gadol-1,k):
n1:=max(k,gadol-1):


gu:=Wilf(n1,Patterns):



S2:={seq([op(max(1,n1-k+1)..n1,w)], w in gu)}:

for w in S2 do
 T4[w]:=0:
od:

for w in gu do
T4[[op(max(1,n1-k+1)..n1,w)]]:=T4[[op(max(1,n1-k+1)..n1,w)]]+1:
od:

S:=[k,S1,op(T1),op(T3),S2,op(T4),n1, WilfSeqBF(Patterns,n1-1)]:
if SchemeSeq(S,n1+1)<>WilfSeqBF(Patterns,n1+1) then
 print(S, `didn't work out`):
 RETURN(FAIL):
else
 RETURN(S):
fi:

end:



#EvalScheme(S,n,L): If a  scheme S=[k,S1,op(T1),op(T3),S2,op(T4),n1]  is supposed to count the number of permutations avoiding the set of patterns Patterns
#then the output is the number of such permutations of length n whose last entries is L. Try:
#S:=FindScheme(7,{[1,2,3]},2,2); EvalScheme(S,8,[6,3]);
EvalScheme:=proc(S,n,L) local k,S1,S2,n1,Avel,i1,a1,pi,V,V1,i,k1,L1:
option remember:
k:=S[1]:
S1:=S[2]:
S2:=S[5]:
n1:=S[7]:

if n<n1 then
  RETURN(FAIL):
fi:

if nops(L)=k and n=n1 then
if not member(L,S2) then
 RETURN(0):
else
  RETURN(S[6][L]):
fi:
fi:


if not member(nops(L),{k,k-1}) then
 RETURN(FAIL):
fi:


if nops(L)=k-1 then
 Avel:={seq(i1, i1=1..n)} minus {op(L)}:
 RETURN(add(EvalScheme(S,n,[a1,op(L)]),a1 in Avel)):
fi:

pi:=redu(L):
if not member(pi,S1) then
  RETURN(0):
fi:

V:=S[4][pi]:

L1:=sort(L):
 V1:=[L1[1]-1,seq(L1[i+1]-L1[i]-1,i=1..k-1),n-L1[k]]:

for i from 1 to nops(V) do
 if V[i]<>1/2 and V1[i]>V[i] then
    RETURN(0):
 fi:
od:



k1:=max(S[3][pi][1]):

L1:=[op(1..k1-1,L),op(k1+1..k,L)]:

L1:=subs({seq(i1=i1-1,i1=L[k1]+1..n)},L1):

EvalScheme(S,n-1,L1):

end:



#EvalSchemeA(S,n): if S is a scheme for counting the number of permutations avoiding Patterns and n is a positive integer it outputs the number of such permuations of length n using the scheme.
EvalSchemeA:=proc(S,n) local k,n1,gu,gu1:
option remember:
k:=S[1]:
n1:=S[7]:

if n<n1 then
 RETURN(FAIL):
fi:

gu:=permute(n,k):

add(EvalScheme(S,n,gu1) ,gu1 in gu):

end:


#WilfSeqBF(Patterns,N): The first N terms, by brute force. Try:
#WilfSeqBF({[1,2,3]},8);
WilfSeqBF:=proc(Patterns,N) local n1: [seq(nops(Wilf(n1,Patterns)),n1=1..N)]:end:


#WilfSeqVS(n,Patterns,k,r,N): first tries to find a scheme for Patterns using  FindSCheme(n,Patterns,k,r) and then uses it to compute the first N terms
#WilfSeqVS(7,{[1,2,3]},2,1,20);
WilfSeqVS:=proc(n,Patterns,k,r,N) local S, n1: 
S:=FindScheme(n,Patterns,k,r,N):
if S=FAIL then
 RETURN(FAIL):
fi:
[op(S[9]),seq(EvalSchemeA(S,n1),n1=S[8]..N)]:
end:




#SchemeSeq(S,N): the first N terms evaluated by the scheme S. Try:
#S:=FindScheme(7,{[1,2,3]},2,1); SchemeSeq(S,20);
SchemeSeq:=proc(S,N) local n1:
[op(S[8]),seq(EvalSchemeA(S,n1),n1=S[7]..N)]:
end: