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


print(`First Written: Feb., 2025: tested for Maple 2020 `):
print(`Version  :  March 11, 2025 `):
print():
print(`This is PPF.txt, A Maple package`):
print(`accompanying Lucy Martinez and Doron Zeilberger's  article: `):
print(` Experimenting with Random Parking Functions`):
print(` http://sites.math.rutgers.edu/~zeilberg/tokhniot/mamarim/mamarimhtml/ppf.html `):
print():
print(`The most current version is available on WWW at:`):
print(` http://sites.math.rutgers.edu/~zeilberg/tokhniot/PPF.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(`-------------------`):


print(`-------------------------------------------`):
print(`For  a list of the Vector Parking functions,`):
print(` type "ezraV();". For specific help type "ezra(procedure_name);" `):
print():
print(`-------------------`):


print(`-------------------------------------------`):
print(`For a list of the Story functions,`):
print(` type "ezraS();". For specific help type "ezra(procedure_name);" `):
print():
print(`-------------------`):



print(`-------------------------------------------`):
print(`For  a list of the generalized Parking functions`):
print(` type "ezraG();". For specific help type "ezra(procedure_name);" `):
print():
print(`-------------------`):


print(`-------------------`):




ezraG:=proc()
if args=NULL then

print(`The Generalized Parking functions procedures are:  BK, ResPP `):
print(``):

else
ezra(args):

fi:


end:


ezraS:=proc()
if args=NULL then

print(`The Story procedures are:  `):
print(``):

else
ezra(args):

fi:


end:


ezra1:=proc()
if args=NULL then

print(`The SUPPORTING procedures are: `):
print(` Ank, CS, CS1, Insert1, NuResPP, ParkingTable, PPgToPPvStupid, PTable, PWilf, Raise1, Tro, TroF, TroV, Wilf, WtEnk `):

else
ezra(args):

fi:


end:



ezraV:=proc()
if args=NULL then

print(`The Vector parking functions procedures are: `):
print(`  CSv, NuPnV, PPv, PnV, RPV `):

else
ezra(args):

fi:


end:


ezra:=proc()
if args=NULL then

print(` PPF.txt: A Maple package for  exploring random parking functions`):

print(`The MAIN procedures are:  CFO, BucketSeqs, NuPna, Park, ParkG, ParkPS, ProbP, PP, PPg, Pna, Tnq, PWilfSeq, WtE, WtPnk  `):
print(` PPgToPPv, MultiFloorSeq, GenerateV`):



elif nargs=1 and args[1]=Ank then
print(`Ank(n,k): [k]^n. Try:`):
print(`Ank(3,3);`):

elif nargs=1 and args[1]=BK then
print(`BK(k,m): The numbers T(k,k*n+s) in the Big buckets are (are not) better!`):
print(` article by Blake and Konheim https://dl.acm.org/doi/pdf/10.1145/322033.322038. Try:`):
print(`seq(BK(2,m),m=1..10);`):

elif nargs=1 and args[1]=CS then
 print(`CS(n): the set of Catalan sequences of length n. Try:`):
 print(`CS(5);`):
print(``):


elif nargs=1 and args[1]=CSv then
 print(`CSv(v): the set of vector Catalan sequences with vector v. Try:`):
 print(`CSv([1,2,3,4,5]);`):
print(``):



elif nargs=1 and args[1]=Insert1 then
print(`Insert1(S,L): Given a set S and a list L such that nops(S)+nops(L)=n, say and S is a subset of {1,...n} first raises all the members of L by 1's`):
print(`and then creates a list of length n where the locations of S are 1's and the other locations follow L. Try:`):
print(`Insert1({2,3},[a,b]);`):


elif nargs=1 and args[1]=NuPna then
print(`NuPna(n,a): the number of  a-parking functions of length n. Try:`):
print(`[seq(NuPna(i,1),i=1..10)];`):

elif nargs=1 and args[1]=NuPnV then
print(`NuPnV(V): the number of  V-parking functions of length n. Try:`):
print(`NuPnV([1,2,3]);`):

elif nargs=1 and args[1]=NuResPP then
print(`NuResPP(n,k,a): the number of  a-parking functions of length n and largest part <=k. Try:`):
print(`[seq(NuResPP(i,i,1),i=1..10)];`):

elif nargs=1 and args[1]=Park then
print(`Park(L): Given  a function L from [1,n] to [1,n] finds the final parking permutation or retunrs FAIL. It also returns the list of where every car winded up Try`):
print(`Park([1,1,1]);`):

elif nargs=1 and args[1]=ParkingTable then
print(`ParkingTable(n): for each permutation the set of parking functions that go with it. Try:`):
print(`ParkingTable(3);`):

elif nargs=1 and args[1]=ParkG then
print(`ParkG(L,C): Given  a function L from [1,n] to [1,n] and parking capacities for each space`):
print(`finds the final parking placement or returss FAIL. It also returns the list of where every car winded up Try`):
print(`ParkG([1,1,1],[1,1,1]);`):

elif nargs=1 and args[1]=ParkPS then
print(`ParkPS(n,t): The weight-enumerator of the permutation according to t^NumberOfParkingFunctionsThatGive(pi)`):
print(`Try:`):
print(`ParkPS(3,t);`):

elif nargs=1 and args[1]=Pna then
print(`Pna(n,a): the set of  a-parking functions of length n. Try:`):
print(`Pna(5,2);`):

elif nargs=1 and args[1]=PP then
 print(`PP(n): the set of parking functions of length n. Try:`):
 print(`PP(5);`):
print(``):

elif nargs=1 and args[1]=PPg then
print(`PPg(C): The set of parking functions with parking capacity C. Try:`):
print(`PPg([1,1,1]);`):


elif nargs=1 and args[1]=PPv then
 print(`PPv(v): the set of vector-parking functions with vector v. Try:`):
 print(`PPv([1,2,3,4,5]);`):
print(``):


elif nargs=1 and args[1]=PnV then
print(`PnV(V): same as PPv(V) but a different approach. Try:`):
print(`PnV([1,2,3]);`):

elif nargs=1 and args[1]=ProbP then
print(`ProbP(P): Given a stochastic matrix P of size n by n (given as list of lists) where P[i][j] is the probability that car i will try to park at location j, outputs`):
print(`the prob. that everyone can park. Try:`):
print(`ProbP([[1/3,1/3,1/6,1/6],[1/3,1/6,1/3,1/6],[1/6,1/3,1/3,1/6],[1/10,1/10,2/5,2/5]]);`):

elif nargs=1 and args[1]=Ptable then
print(`Ptable(pi) : the parking table of the permutation pi. Try:`):
print(`Ptable([4,1,2,3]);`):

elif nargs=1 and args[1]=PWilf then
print(`PWilf(Patterns,n): The number of parking functions that yield permutations avoiding the patterns in the set of patterns Patterns. Try:`):
print(`PWilf({[1,2,3],[1,3,2]},n);`):

elif nargs=1 and args[1]=PWilfSeq then
print(`PWilfSeq(N,Patterns): inputs a positive integer N, and a set of patterns Patterns, and outputs the list of length N whose n-th entry is`):
print(`the number of parking functions of length n that yield permutations avoiding the patterns in the set of patterns Patterns. Try:`):
print(`PWilfSeq({[1,2,3],[1,3,2]},6);`):

elif nargs=1 and args[1]=Raise1 then
print(`Raise1(f,a,n,m): Given a polynomial f in a[i,j] where 1<=i<=n and 1<=j<=m, replaces a[i,j] by a[i,j+1]. Try:`):
print(`Raise1(a[1,3]*a[2,4],a,2,7);`):

elif nargs=1 and args[1]=ResPP then
print(`ResPP(n,k) he set of restricted parking unctions of length n with largest entry k.`):
print(`ResPP(n,,n) is the same as PP(n). Try:`):
print(`ResPP(5,2);`):


elif nargs=1 and args[1]=Tnq then
print(`Tnq(n,q): the trouble polynomial. Try:`):
print(`Tnq(5,q);`):

elif nargs=1 and args[1]=Tro then
print(`Tro(L): the total trouble of L for a usual parking finction. Try`):
print(`Tro([1,1,1]);`):


elif nargs=1 and args[1]=TroF then
print(`TroF(L): Same as Tro(L) but faster. Try`):
print(`TroF([1,1,1]);`):

elif nargs=1 and args[1]=TroV then
print(`TroV(L): the total trouble vector of parking function L. Try:`):
print(`TroV([2,2,1,2]);`):



elif nargs=1 and args[1]=Wilf then
print(`Wilf(n,Patterns): all permutations of length n`):
print(`that avoid the patterns in the set of patterns`):
print(`Patterns. Try: `):
print(`Wilf(5,{[1,2,3],[1,3,2]});`):
 
elif nargs=1 and args[1]=WtE then
print(`WtE(n,a): the weight-enumertor of parking function done directly. Try:`):
print(`WtE(5,a);`):

elif nargs=1 and args[1]=WtEnk then
print(`WtEnk(n,k,a): the weight-enumertor of k-parking function done directly`):

elif nargs=1 and args[1]=WtPnk then
print(`WtPnk(n,k,a): Same as WtEnk(n,k,a) but cleverly. Try:`):
print(`WtPnk(5,1,a);`):

elif nargs=1 and args[1]=Wt1 then
print(`Wt1(L,a): the weight of the parking function L w.r.t. a[i,p[i]] `):



elif nargs=1 and args[1]=PPgToPPv then
print(`PPgToPPv(L): inputs a list of pos. integers that denote the number of`):
print(`floors in each parking place and outputs, if possible,`):
print(`a vector V such that it equals PPv(V). Try:`):
print(`PPgToPPv([2,2,2]);`):

elif nargs=1 and args[1]=PPgToPPvStupid then
print(`PPgToPPvStupid(L): inputs a list of pos. integers that denote the number of`):
print(`floors in each parking place and outputs, if possible,`):
print(`a vector V such that it equals PPv(V). Try:`):
print(`PPgToPPvStupid([2,2,2]);`):


elif nargs=1 and args[1]=CFO then
print(`CFO(outcome): Given a permutation outcome, counts the number of`):
print(`parking functions that park in the order of the permutation outcome using`):
print(`the counting through permutations method. Try:`):
print(`CFO([1,2,3,4]);`):



elif nargs=1 and args[1]=BucketSeqs then
print(`BucketSeqs(k,N): Inputs a positive integer k and`):
print(`a positive integer N, and outputs a list of k lists`):
print(`whose i-th list is the sequence T(k,kn+i-1) of the great paper by I.A. Blake and A.G. Konheim J. ACM, 24 (1977), 591-606.`):
print(`For example to get the first 10 terms of OEIS sequence https://oeis.org/A006698 `):
print(`and https://oeis.org/A068087 , Try:`):
print(`BucketSeqs(2,10);`):
print(`Yet another example, to get`):
print(`https://oeis.org/A006699, https://oeis.org/A006700, and https://oeis.org/A208148, try:`):
print(`BucketSeqs(3,10);`):

elif nargs=1 and args[1]=MultiFloorSeq then
print(`MultiFloorSeq(k,N): inputs pos. integers k and N and outputs a list L of length N`):
print(`such that L is exactly the list [seq(nops(PPg([k$n]), n=1..N)] which is`):
print(`the set of parking functions with parking capacity C. Try:`):
print(`MultiFloorSeq(1,12);`):


elif nargs=1 and args[1]=GenerateV then
print(`GenerateV(n,k): Generates a vector V of length n for which PPv(V)=PPg([k,k,k,...,k])`):
print(`where k is repeated n times. Try:`):
print(`GenerateV(6,2);`):

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

fi:


end:


with(Statistics):
##start simulation
HelpSimu:=proc()

if args=NULL then

print(`The  SIMULATION procedures are:  CheckPF , EstimatePr, EstProbP, IsP,  Roullete, RPM, RPna, RPV, RV `):


elif nargs=1 and args[1]=EstProbP then
print(`EstProbP(P,K): Given a stocahstic matrix P of size n by n (given as list of lists) where P[i][j] is the probability that car i will try to park at location j, estimates ProbP(P)`):
print(`by averaging over K random parking functions. Try:`):
print(`P:=RPM(4,20): EstProbP(P,100);`):

elif nargs=1 and args[1]=IsP then
print(`IsP(P): Inputs an n by n probability matrix P=[[p11,p12,...,p1n],[p21,p22,...,p2n],...,[pn1,pn2,...,pnn]] `):
print(`where pi1+pi2+...+pin=1 for all 1<=i<=n, and checks that`):
print(`all probabilities add up to 1 and that each list is size n, returns true or false`):

elif nargs=1 and args[1]=RPM then
print(`RPM(n,K): random n by n prob. matrix. Try:`):
print(`RPM(20,100);`):

elif nargs=1 and args[1]=RV then
print(`RV(P): inputs a probability matrix such that P[i][j] is the probability that`):
print(`car i gets assigned spot j and outputs the preference vector of length n. Try:`):
print(`RV([[1/3,2/3],[1/4,3/4]]); `):


elif nargs=1 and args[1]=CheckPF then
print(`CheckPF(a): inputs a preference vector a of length n and checks whether a is`):
print(`a parking function via the inequality condition. Try:`):
print(`CheckPF([1,1,1]); `):



elif nargs=1 and args[1]=EstimatePr then
print(`EstimatePr(P,K): inputs an integer n, a probability matrix P and an integer K,`):
print(`runs RV(P) "K" times, and outputs the ratio of times where the random vector was a parking function. Try:`):
print(`EstimatePr([[1/3,1/3,1/3],[1/3,1/3,1/3],[1/3,1/3,1/3]],1000); `):

elif nargs=1 and args[1]=Roullete then
print(`Roullete(L): Given a list of positive integers L outputs i with prob, L[i]/add(L). Try:`):
print(`Roullete ([1,2,3]);`):


elif nargs=1 and args[1]=RPna then
print(`RPna(n,a): a random a-parking function of length n. Try: `):
print(`RPna(10,1);`):


elif nargs=1 and args[1]=RPV then
print(`RPV(V): a random V-parking function with vector V. Try`):
print(`RPV([1,2,3,4,5]);`):

elif nargs=1 and args[1]=SpecialPP then
print(` SpecialPP(n,x,K): inputs a number between 0 and 1 and`):
print(`  EstimatePr([[x,(1-x)/(n-1)]$n],K) `):
print(`  Then experiment with a large K, say K=10000`):
print(`   and n=40 and plot it for x from 0 to 1 with`):
print(`   x=0.01*i, i=1..100, then plot it.`):

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

end:


###############################
#IsP(P): Inputs an n by n probability matrix P=[[p11,p12,...,p1n],[p21,p22,...,p2n],...,[pn1,pn2,...,pnn]]
# where pi1+pi2+...+pin=1 for all 1<=i<=n,
# and checks that all probabilities add up to 1 and that each list is size n, returns
# true or false
IsP:=proc(P) local s,i,n:
n:=nops(P):

for s in P do
 if not (add(i, i in s)=1) then
  print(`Sum of all the probabilities should add up to 1`):
  print(`Check the following table`, s):
  return false:
 fi:

 if not nops(s)=n then
  print(`The size of each list should be`, n):
  return false:
 fi:
od:
true:
end:


#RV(P): inputs a probability matrix such that P[i][j] is the probability that
# car i gets assigned spot j and outputs the preference vector of length n
RV:=proc(P) local X,X1,s,a:
IsP(P):
a:=[]:
for s in P do
 X:=RandomVariable(ProbabilityTable(s)):
 #X1 returns 1,2,..., or n
 X1:=Sample(X,1):
 a:=[op(a),round(X1[1])]:
od:
a:
end:


#CheckPF(a): inputs a preference vector a of length n and checks whether a is
# a parking function via the inequality condition
CheckPF:=proc(a) local n,i,newa:
n:=nops(a):
newa:=sort(a):
for i from 1 to n do
 if not newa[i]<=i then
  return false:
 fi:
od:
true:
end:

# EstimatePr(P,K): inputs an integer n, a probability matrix P and an integer K,
#  runs RV(P) "K" times, and outputs the ratio of times where
#  the random vector was a parking function
EstimatePr:=proc(P,K) local co,a,i:
co:=0:
for i from 1 to K do
 #generate a random vector a
 a:=RV(P):
 if CheckPF(a) then
  co:=co+1:
 fi:
od:
evalf(co/K):
end:


# SpecialPP(n,x,K): inputs a number between 0 and 1 and
#  EstimatePr([[x,(1-x)/(n-1)]$n],K)
#  Then experiment with a large K, say K=10000
#   and n=40 and plot it for x from 0 to 1 with
#   x=0.01*i, i=1..100, then plot it.
SpecialPP:=proc(n,x,K) local P,i:
P:=[[x,seq((1-x)/(n-1),i=1..n-1)]$n]:
EstimatePr(P,K):
end:
##End simulation


with(combinat):

#CS1(n,k):
CS1:=proc(n,k) local gu,lu,lu1,k1:
option remember:
if n=1 then
if k=1 then
 RETURN({[1]}):
else
RETURN({}):
fi:
fi:

gu:={}:

for k1 from 1 to k do
lu:=CS1(n-1,min(k1,n-1)):
gu:=gu union {seq([op(lu1),k],lu1 in lu)}:
od:
gu:
end:

#CS(n): the set of Catalan sequences of length n. 
CS:=proc(n) local k:
option remember:
{seq(op(CS1(n,k)),k=1..n)}:
end:

#PP(n): the set of parking functions of length n
PP:=proc(n) local gu,gu1:
gu:=CS(n):
{seq(op(permute(gu1)),gu1 in gu)}:
end:


#Wt1(L,a): the weight of the parking function L
Wt1:=proc(L,a) local i:
mul(a[i,L[i]],i=1..nops(L)):
end:

#WtE(n,a): the weight-enumerator of parking function done directly
WtE:=proc(n,a) local gu,gu1:
gu:=PP(n):
add(Wt1(gu1,a), gu1 in gu):
end:

#WtEnk(n,k,a): the weight-enumertor of k-parking function done directly
WtEnk:=proc(n,k,a) local gu,gu1:
gu:=Pna(n,k):
add(Wt1(gu1,a), gu1 in gu):
end:

#ProbP(P): Given a stochastic matrix P of size n by n (given as list of lists) where P[i][j] is the probability that car i will try to park at location j, outputs
#the prob. that everyone can park
ProbP:=proc(P) local n,i,gu,a,j:
n:=nops(P):
if {seq(nops(P[i]),i=1..n)}<>{n} and {seq(convert(P[i],`+`),i=1..n)}<>{1} then
 RETURN(FAIL):
fi:
gu:=WtE(n,a):
subs({seq(seq(a[i,j]=P[i][j],j=1..n),i=1..n)},gu):
end:


#EstProbP(P,K): Given a stocahstic matrix P of size n by n (given as list of lists) where P[i][j] is the probability that car i will try to park at location j, estimates ProbP(P)
#by averaging over K random parking functions. Try:
#P:=RPM(4,20): EstProbP(P,100);
EstProbP:=proc(P,K) local n,i,gu,L,i1:
n:=nops(P):
if {seq(nops(P[i]),i=1..n)}<>{n} and {seq(convert(P[i],`+`),i=1..n)}<>{1} then
 RETURN(FAIL):
fi:
gu:=0:

for i from 1 to K do
L:=RPna(n,1):
gu:=gu+mul(P[i1][L[i1]],i1=1..n):
od:
evalf(gu/K*(n+1)^(n-1)):
end:




#Insert1(S,L): Given a set S and a list L such that nops(S)+nops(L)=n, say and S is a subset of {1,...n} first raises all the members of L by 1's
#and then creates a list of length n where the locations of S are 1's and the other locations follow L. Try:
#Insert1({2,3},[a,b]);
Insert1:=proc(S,L) local n,i,M,L1:
n:=nops(S)+nops(L):
if S minus {seq(i,i=1..n)}<>{} then
 RETURN(FAIL):
fi:
L1:=L+[1$nops(L)]:
M:=[]:
for i from 1 to n do
 if member(i,S) then
  M:=[op(M),1]:
 else
 M:=[op(M),L1[1]]:
 L1:=[op(2..nops(L1),L1)]:
fi:
od:
M:
end:


#Pna(n,a): the set of  a-parking functions of length n
Pna:=proc(n,a) local gu,lu,lu1,r,mu,mu1:
option remember:

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

gu:={}:

lu:=powerset(n):

for lu1 in lu do
r:=nops(lu1):
mu:=Pna(n-r,a+r-1):
gu:=gu union {seq(Insert1(lu1,mu1),mu1 in mu)}:
od:
gu:
end:


#Raise1(f,a,n,m): Given a polynomial f in a[i,j] where 1<=i<=n and 1<=j<=m, replaces a[i,j] by a[i,j+1]. Try:
#Raise1(a[1,3]*a[2,4],a,2,7);
Raise1:=proc(f,a,n,m) local i,j:
subs({seq(seq(a[i,j]=a[i,j+1],j=1..m),i=1..n)},f):
end:




#WtPnk(n,k,a): Same as WtEnk(n,k,a) but cleverly. Try:
#WtPnk(5,1,a);
WtPnk:=proc(n,k,a) local gu,lu,lu1,r,mu,lu1c,i1,j1:
option remember:

if k=0 then
 RETURN(0):
fi:
if n=0 then
RETURN(1):
fi:

gu:=0:

lu:=powerset(n):

for lu1 in lu do
r:=nops(lu1):

lu1c:=sort(convert({seq(i1,i1=1..n)} minus lu1,list)):

mu:=WtPnk(n-r,k+r-1,a):
mu:=Raise1(subs({seq(seq(a[i1,j1]=a[lu1c[i1],j1],i1=1..n-r),j1=1..n+k)},mu),a,n,n+k):

gu:=expand(gu+ mul(a[i1,1], i1 in lu1)*mu):
od:
gu:
end:

#Park(L): Given  a function L from [1,n] to [1,n] finds the final parking permutation or returns FAIL. It also returns the list of where every car winded up Try
#Park([1,1,1]);
Park:=proc(L) local n,L1,i,j,a,L2:
n:=nops(L):

if {op(L)} minus {seq(i,i=1..n)}<>{} then
 ERROR(`bad input`):
fi:

L1:=[0$n]:

for i from 1 to n do
 a:=L[i]:

 if L1[a]=0 then
   L1:=[op(1..a-1,L1),i,op(a+1..n,L1)]:
 L2[i]:=a:
 else
 for j from a+1 to n while L1[j]<>0 do od:
  if j=n+1 then
    RETURN(FAIL):
 else
   L1:=[op(1..j-1,L1),i,op(j+1..n,L1)]:
 L2[i]:=j:
 fi:
fi:

od:

[L1,[seq(L2[j],j=1..n)]]:

end:


#TroV(L): the total trouble vector of parking function L. Try:
#TroV([2,2,1,2]);
TroV:=proc(L) local n,L1,i:
if not CheckPF(L) then
 RETURN(FAIL):
fi:

n:=nops(L):
L1:=Park(L)[2]:
[seq(L1[i]-L[i],i=1..n)]:
end:

#Tro(L): the total trouble of L
Tro:=proc(L) :
convert(TroV(L),`+`):
end:



#TroF(L): the total trouble of L
TroF:=proc(L) :
binomial(nops(L)+1,2)-convert(L,`+`):
end:


#Tnq(n,q): the trouble polynomial. Try:
#Tnq(5,q);
Tnq:=proc(n,q) local gu,gu1:
gu:=PP(n):
add(q^Tro(gu1),gu1 in gu):
end:





#ParkG(L,C): Given  a function L from [1,n] to [1,n] and parking capacities for each space
#finds the final parking placement or retunrs FAIL. It also returns the list of where every car winded up Try
#ParkG([1,1,1],[1,1,1]);
ParkG:=proc(L,C) local n,L1,i,j,a,L2:
n:=nops(L):

if not (type(C,list) and nops(C)=n) then
 ERROR(`bad input`):
fi:

if {op(L)} minus {seq(i,i=1..n)}<>{} then
 ERROR(`bad input`):
fi:

L1:=[[]$n]:

for i from 1 to n do
 a:=L[i]:

 if nops(L1[a])<C[a] then
   L1:=[op(1..a-1,L1),[op(L1[a]),i],op(a+1..n,L1)]:
 L2[i]:=a:
 else
 for j from a+1 to n while nops(L1[j])=C[j] do od:
  if j=n+1 then
    RETURN(FAIL):
 else
   L1:=[op(1..j-1,L1),[op(L1[j]), i],op(j+1..n,L1)]:
 L2[i]:=j:
 fi:
fi:

od:

[L1,[seq(L2[j],j=1..n)]]:

end:



#Ank(n,k): [k]^n
Ank:=proc(n,k) local gu,i,gu1:
option remember:
if n=0 then
 RETURN({[]}):
fi:
gu:=Ank(n-1,k):
{seq(seq([op(gu1),i],i=1..k),gu1 in gu)}:
end:


#PPg(C): The set of parking functions with parking capacity C. Try:
#PPg([1,1,1]);
PPg:=proc(C) local n,mu,mu1,gu:
n:=nops(C):
mu:=Ank(n,n):
gu:={}:
for mu1 in mu do
 if ParkG(mu1,C)<>FAIL then
  gu:=gu union {mu1}:
 fi:
od:
gu:
end:



#RP1(n,K):  a random prob. vector of length n using K
RP1:=proc(n,K) local ra,L,i:
ra:=rand(1..K):
L:=[seq(ra(),i=1..n)]:
L/add(L):
end:


#RPM(n,K): random n by n prob. matrix
RPM:=proc(n,K) local i:
[seq(RP1(n,K),i=1..n)]:
end:


#Roullete(L): Given a list of positive integers L outputs i with prob, L[i]/add(L). Try:
#Roullete ([1,2,3]);
Roullete:=proc(L) local i,j,r:

if not (type(L,list) and {seq(type(L[i],integer),i=1..nops(L))}={true} and {seq(evalb(L[i]>=0),i=1..nops(L))}={true}) then
 RETURN(FAIL):
fi:

r:=rand(1..add(L))():
for i from 1 while r>add(L[j],j=1..i) do 

od:
i:


end:


#RPna(n,a): a random a-parking function of length n. Try: 
#RPna(10,1);
RPna:=proc(n,a) local r,L,gu,mu,S,i:
if n=0 then
 RETURN([]):
fi:

L:=[seq(binomial(n,r)*(a+r-1)*(a+n-1)^(n-r-1),r=0..n)]:

r:=Roullete(L)-1:


S:=randcomb(n,r):


gu:=RPna(n-r,a+r-1)+[1$(n-r)]:

mu:=[]:

for i from 1 to n do
if member(i,S) then
 mu:=[op(mu),1]:
else
 mu:=[op(mu),gu[1]]:
 gu:=[op(2..nops(gu),gu)]:
fi:
od:

mu:
end:



#CSv1(v,k): the increasing sequences a such that a[i]<=v[i] and a[nops(v)]=k
CSv1:=proc(v,k) local n, gu,lu,k1,lu1:
option remember:
n:=nops(v):
if n=1 then
 if k<=v[1] then
 RETURN({[k]}):
else
 RETURN({}):
fi:
fi:

if  not (k>=1 and k<=v[n]) then
 RETURN({}):
fi:

gu:={}:

for k1 from 1 to min(k,v[n-1]) do
lu:=CSv1([op(1..n-1,v)],k1):
gu:=gu union {seq([op(lu1),k],lu1 in lu)}:
od:
gu:
end:

CSv:=proc(v) local n,k:
option remember:
n:=nops(v):
{seq(op(CSv1(v,k)),k=1..v[n])}:
end:




#PPv(v): the set of vector parking functions with vector v
PPv:=proc(v) local gu,gu1:
gu:=CSv(v):
{seq(op(permute(gu1)),gu1 in gu)}:
end:






#PnV(n,V): the set of  V-parking functions of length n. Try:
#PnV(3,[1,2,3]);
PnV:=proc(V) local n,gu,lu,lu1,r,mu,mu1:
option remember:

gu:={}:

n:=nops(V):

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

if member(0,V) then
 RETURN({}):
fi:


lu:=powerset(n):

for lu1 in lu do
r:=nops(lu1):
mu:=PnV([op(r+1..n,V)]-[1$(n-r)]):
gu:=gu union {seq(Insert1(lu1,mu1),mu1 in mu)}:
od:
gu:
end:



#NuPnVold(V): the number of  V-parking functions of length n. Try:
#NuPnVold([1,2,3]);
NuPnVold:=proc(V) local n,gu,lu,lu1,r:
option remember:

gu:=0:

n:=nops(V):

if n=0 then
RETURN(1):
fi:

if member(0,V) then
 RETURN(0):
fi:


lu:=powerset(n):

for lu1 in lu do
r:=nops(lu1):
gu:=gu+NuPnV([op(r+1..n,V)]-[1$(n-r)]):
od:
gu:
end:





#ResPP(n,k) the set of restricted parking unctions of length n with largest entry k.
#ResPP(n,n) is the same as PP(n). Try:
#ResPP(5,2);
ResPP:=proc(n,k) local gu,gu1,k1:
gu:={seq(op(CS1(n,k1)),k1=1..min(k,n))}:
{seq(op(permute(gu1)),gu1 in gu)}:
end:



#NuPna(n,a): the number of  a-parking functions of length n. Try:
#[seq(NuPna(i,1),i=1..10)];
NuPna:=proc(n,a) local r:
option remember:

if a=0 then
 RETURN(0):
fi:
if n=0 then
RETURN(1):
else(add(binomial(n,r)*NuPna(n-r,a+r-1),r=0..n)):
fi:

end:


#NuResPP(n,k,a): the number of  a-parking functions of length n and largest part <=k. Try:
#[seq(NuResPP(i,i,1),i=1..10)];
NuResPP:=proc(n,k,a) local r:
option remember:


if n=0 then
RETURN(1):
fi:

if a=0 or k=0 then
 RETURN(0):
fi:

add(binomial(n,r)*NuResPP(n-r,k-1,a+r-1),r=0..n):

end:

#PPgToPPvStupid(L): inputs a list of pos. integers that denote the number of
#floors in each parking place and outputs, if possible,
#a vector V such that it equals PPv(V)
PPgToPPvStupid:=proc(L) local S,S1,S2,s,s1,V,i:
S:=PPg(L):
S1:={seq(sort(s),s in S)}:
S2:={}:

#Check that if you permute all of S1 and take the union you get back S
for s in S1 do
 s1:=convert(permute(s),set):
 S2:=S2 union s1:
od:

if S minus S2 <> {} then
 return(FAIL):
fi:

#For i from 1 to n, finds the minimum of s1[i] over all s1 in S1
V:=[seq(max(seq(s1[i], s1 in S1)), i = 1 .. nops(s1[1]))]:
if S=PPv(V) then
 return V:
else
 return(FAIL):
fi:
end:


#PPgToPPv(L): Like PPgToPPvStupid(L), but done cleverly
PPgToPPv:=proc(L) local n,V,L1,n1:
L1:=L:
V:=[]:
n:=nops(L):
n1:=nops(L):
while nops(V)<n and L1<>[] do
V:=[  n1$L1[-1]  ,op(V)]:
n1:=n1-1:
L1:=[op(1..nops(L1)-1,L1)]:
od:
if nops(V)>n then
V:=[op(nops(V)-n+1..nops(V),V)]:
fi:
V:
end:



#NuPnV(V): the number of  V-parking functions of length n. Try:
#NuPnV([1,2,3]);
NuPnV:=proc(V) local n,gu,r:
option remember:

gu:=0:

n:=nops(V):

if n=0 then
RETURN(1):
fi:

if member(0,V) then
 RETURN(0):
fi:


for r from 0 to n do
gu:=gu+binomial(n,r)*NuPnV([op(r+1..n,V)]-[1$(n-r)]):
od:
gu:
end:




#RPV(V): a random V-parking function with vector V. Try
#RPV([1,2,3,4,5]);
RPV:=proc(V) local n,r,L,gu,mu,S,i:

n:=nops(V):
if n=0 then
 RETURN([]):
fi:

L:=[seq(binomial(n,r)*NuPnV([op(r+1..n,V)]-[1$(n-r)]),r=0..n)]:
r:=Roullete(L)-1:


S:=randcomb(n,r):


gu:=RPV([op(r+1..n,V)]-[1$(n-r)])+[1$(n-r)]:

mu:=[]:

for i from 1 to n do
if member(i,S) then
 mu:=[op(mu),1]:
else
 mu:=[op(mu),gu[1]]:
 gu:=[op(2..nops(gu),gu)]:
fi:
od:

mu:
end:


#ParkingTable(n): for each permutation the set of parking functions that go with it
ParkingTable:=proc(n) local gu,mu,T,i,mu1:
gu:=permute([seq(i,i=1..n)]):
mu:=PP(n):

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

for mu1 in mu do
T[Park(mu1)[1]]:=T[Park(mu1)[1]] union {mu1}:
od:

[seq([gu[i],T[gu[i]]],i=1..nops(gu))];

end:

#ParkPS(n,t): The weight-enumerator of the permutation according to t^NumberOfParkingFunctionsThatGive(pi)
#Try:
#ParkPS(3,t);
ParkPS:=proc(n,t) local gu,i:
gu:=ParkingTable(n):
add(t^nops(gu[i][2]),i=1..nops(gu)):
end:


#CFOold(outcome): Given an permutation outcome, counts the number of parking functions that
#park in the order of the permutation outcome using the counting through permutations method
#Try: CFOold([3,2,1]);
CFOold:=proc(outcome) local n, i, j, Li, product:
n:=nops(outcome):
if {op(outcome)}<>{seq(i,i=1..n)} then
 print(`The given input must be a permutation`):
 return(FAIL):
fi:

product:=1: 
for i from 1 to n do
 Li :=0:  #start with a length of 1 for L(i, perm)
 #compute L(i,outcome) by looking backwards
 for j from i by -1 to 1 do
  if outcome[j]<=outcome[i] then
   Li := Li + 1:
  else
   break: #stop counting when the condition is violated
  fi:
 od:
 product:=product*Li:
od:
product:
end:


#BucketSeqs(k,N): Inputs a positive integer k and
BucketSeqs:=proc(k,N) local i,j,n:
{seq([seq(NuPnV([seq(i$k,i=1..n-1),n$j]),n=1..N)],j=0..k-1)}
end:



##############################Start from Wilf


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


#good1(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 0 if this is
#the case, otherwise 1 (0 means bad)

good1:=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(1):
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(0):
fi:
od:

1:

end:



#good(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 0 if this is
#the case, otherwise 1 (0 means bad)

good:=proc(pi,SetPatterns)
local i:
if member(pi,SetPatterns) then
RETURN(0):
fi:

for i from 1 to nops(SetPatterns) do
if good1(pi,op(i,SetPatterns))=0 then
 RETURN(0):
fi:
od:
1:
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 good(muamad,Patterns)=1 then
   gu:=gu union {muamad}:
  fi:
fi:
od:
od:
 
gu:
end:

#PWilf(Patterns,n): The number of parking functions that yield permutations avoiding the patterns in the set of patterns Patterns. Try:
#PWilf({[1,2,3],[1,3,2]},n);
PWilf:= proc(Patterns,n) local gu, gu1:
gu:=Wilf(n, Patterns):
add(CFO(gu1), gu1 in gu):
end:


#PWilfSeq(N,Patterns): inputs a positive integer N, and a set of patterns Patterns, and outputs the list of length N whose n-th entry is
#the number of parking functions of length n that yield permutations avoiding the patterns in the set of patterns Patterns. Try:
#PWilfSeq({[1,2,3],[1,3,2]},6);
PWilfSeq:=proc(Patterns, N) local n:
seq(PWilf(Patterns,n), n = 1 .. N):
end:

###end from Wilf
#GenerateV(n,k): Generates a vector V of length n for which PPv(V)=PPg([k,k,k,...,k])
# where k is repeated n times
GenerateV:=proc(n,k) local v,j,i;
v:=[]:
j:=n:
#iterate the loop until we have a vector of length n
while nops(v)<n do
 #add k duplicates of n, n-1,n-2, and so on from the end of the vector
 # to the beginning until we have
 # a vector of lenght n. Take the minimum of k and n-nops(v) when we get to the the beginning
 v:=[seq(j,i=1..min(k, n-nops(v))),op(v)]:
 j:=j-1:
od:
v:
end:

#MultiFloorSeq(k,N): inputs pos. integers k and N and outputs a list L of length N
# such that L is exactly the list [seq(nops(PPg([k$n]), n=1..N)] which is
# the set of parking functions with parking capacity C
MultiFloorSeq:=proc(k,N) local n:
[seq(NuPnV(GenerateV(n,k)),n=1..N)]:
end:


#BK(k,m): The numbers T(k,k*n+s) in the Big buckets are (are not) better!
# article by Blake and Konheim https://dl.acm.org/doi/pdf/10.1145/322033.322038
BK:=proc(k, m) local n, i,lu:
option remember:
if m <= k then
 return 1:
else
 n:=iquo(m, k);
 lu:=add(binomial(m - 1,k*i - 1)*i*BK(k,k*i-1)*BK(k,m-k*i),i=1..n):
 if n*k=m then
  RETURN(lu):
 else
  RETURN(lu+(n + 1)*BK(k, m - 1)):
 fi:
fi:
end:


#CFOrec(pi): recursive version of CFO(pi)
CFOrec:=proc(pi) local i,n: 
option remember:
n:=nops(pi):
if pi=[1] then  
RETURN(1): 
fi:

for i from 0 to n-1 while pi[n-i]<=pi[n] do od:
i*CFO(redu([op(1..n-1,pi)])):
end:



#CFO(outcome): Given a permutation outcome, counts the number of
#parking functions that park in the order of the permutation outcome using
#the counting through permutations method
CFO:=proc(pi) local i,j,n,p: 
n:=nops(pi):
p:=1:
for j from 1 to n do
for i from 0 to j-1 while pi[j-i]<=pi[j] do od:
p:=i*p:
od:
end:




#Ptable(pi) : the parking table of the permutation pi. Try:
#Ptable([4,1,2,3]);
Ptable:=proc(pi) local i,j,n,L: 
n:=nops(pi):
L:=[]:
for j from 1 to n do
for i from 0 to j-1 while pi[j-i]<=pi[j] do od:
L:=[op(L),i]:
od:
L:
end: