###ParkingAndTrees.txt###


##########################################################################
## ParkingAndTrees.txt Save this file as ParkingAndTrees.txt            ##
## to use it, stay in the same directory, get into Maple                ##
## (by typing: maple <Enter> ) and then type:                           ##
## read ParkingAndTrees.txt <Enter>                                     ##
## Then follow the instructions given there                             ##
##                                                                      ##
## Written by Yukun Yao and Doron Zeilberger, Rutgers University ,      ## 
## yao@math.rutgers.edu and DoronZeil@gmail.com                         ##
##########################################################################

with(combinat):
with(plots):

print(`First Written: April 2018`):
print(`Version 2018:`):
print():
print(`This is ParkingAndTrees.txt, one of the  Maple packages`):
print(`accompanying Yukun Yao’s and Doron Zeilberger’s  article: `):
print(`An Experimental Mathematics Approach to the Sum of Parking Functions `):

print():
print(`The most current version is available on WWW at:`):
print(` https://sites.math.rutgers.edu/~yao/ParkingAndTees.txt .`):
print(`Please report all bugs to: yao at math dot rutgers dot edu .`):
print():
print(`For an overview of the MAIN functions type “Help()”`):
print(`For a list of the supporting functions type: ezra1();`):
print(`For a list of the Checking procedures type ezraC();, for help with`):
print(`For general help, and a list of the MAIN functions,`):
print(` type "ezra();". For specific help type "ezra(procedure_name);" `):
print():
 
Help:=proc(): 
print(` Cn(n): the set of weakly increasing parking functions of length n.`):
print(` Pn(n): the set of parking functions of length n. `):
print(` cn(n): the number of weakly increasing parking functions of length n obtained by recurrence relations.`): 
print(` pn(n): the number of parking functions of length n obtained by recurrence relations.`):
print(` Cna(n,a): the set of weakly increasing a-parking functions of length n, i.e. the ith element <=a+i-1.`):
print(` Pna(n,a): the set of a-parking functions of length n.`):
print(` PnaNew(n,a): the set of a-parking functions. New format: [a, P] where P is the a-parking function. `):
print(` cna(n,a): the number of weakly increasing a-parking functions of length n obtained by recurrence relations.`): 
print(` pna(n,a): the number of a-parking functions of length n obtained by recurrence relations.`): 
print(` Tna(n,a): the set of rooted labelled forests on (a+n) vertices with a components, where the roots must be labelled from 1 to a.`):
print(` PtoT(L,a): map an a-parking function L (as a list) to a forest in Tna(n,a) non-recursively. `):
print(` TtoP(T): map a forest T in Tna(n,a) to an a-parking function non-recursively. `):
print( `TnaGeneral(n,a): the set of forests of rooted labelled trees on (a+n) vertices with a components, where there is no restriction on which vertices can be the roots.`):
print(` Inv(T): input a forest in TnaGeneral(n,a) (if connected then it is a tree, i.e. a=1), output the inversion number of the forest.`):
print(` PairDist(T): input a tree, output the total pair distances over all pairs of its vertices. `):
print(` PnzT(z,n): the sum of z^Inv(T) over all labelled rooted trees on [n].`):
print(` PnzF(z,n): the sum of z^Inv(F) over all labelled rooted forests on [n].`):
print(` Pnax(n,a,x): the sum of x^(sum of all entries in the parking function) over the set of a-parking functions of length n by recurrence relation.`):
print(` PnPD(z,n): the sum of z^PairDist(T) over all labelled rooted trees on [n].`):
print(` SumOfHeight(T): input a tree or forest, output the sum of heights of all vertices.`):
print(` TotalHeightT(x,n): the sum of x^SumOfHeight(T) over all labelled rooted trees on [n].`):
print(` TotalHeightF(x,n): the sum of x^SumOfHeight(T) over all labelled rooted forests on [n].`):
print(` ImagePtoT(n,a): output the image of all a-parking functions of length n under the map PtoT.`):
print(` ImagePtoTall(n,a): output the image of all a-parking functions of length n under the map PtoT and the pre-image together.`):
print(` ImagePtoTboth(n): output the image of all 1-parking functions of length n under the map PtoT, and PtoT0, and the pre-image together.`):
print(` ImageTtoP(n,a): output the image of all labelled rooted forests on [a+n] with a components where [a] are the a roots under the map TtoP.`):
print(` ImageTtoPboth(n,a): output the image of all labelled rooted forests on [a+n] with a components where [a] are the a roots under the map TtoP and the pre-image together.`):
print(` TnEB(n): the set of labelled rooted trees on [n] by considering which vertex is the root and using Tna1(Emps, Bos).`):
print(` FnEB(n): the set of labelled rooted forests on [n] by considering which vertices are the roots and using Tna1(Emps, Bos).`):
print(` CanF(T,r): the canonical form of a rooted tree T with root r.`):
print(` DrawRT(T,r): draws the rooted tree T with root r.`):
print(` Milon(n): all the images of the 1-parking functions of length n as rooted trees with 0 the root.`):
print(` PtoT0(P): inputs a 1-parking function and outputs the rooted tree with root 0 and non-roots labelled {1, ..., nops(P)}. `):
print(` RPF(n,a): a random a-parking function of length n. `):
print(` RRT(n): a random labelled rooted tree with root 0 and non-roots labelled {1, ..., n}. `):
print(` tn(n,a): the number of forests of rooted trees  where each component has the form [root, SetOfEdges] n employees and a bosses, i.e. the number of ordered labelled forests. `):
print(` TnaEB(Emps,Bos): The set of sets of rooted trees where the set of  roots is Bos and the set of non-root vertices is Emps. `):
print(` Tn0(n): The set of rooted trees on {0,...,n} rooted at 0.`):
print(` PtoTEB(P,Emps,Bos): inputs an a:=nops(Bos)-parking function and outputs the corresponding forest of rooted trees with a components, with a prescribed set of Employers and Bosses, Emps and Bos resepctively. This map is defined recursively.`):

end:

ezraC:=proc()
if args=NULL then
 print(` The Checking procedures are: CheckfP, CheckgP, CheckPtoTEB, CheckTtoP, CheckPtoT, CheckPinverse, CheckTinverse `):
 print(``):
else
ezra(args):
fi:
end:

ezra1:=proc()
if args=NULL then
 print(` The supporting procedures are: Cna1, Conv1, Dorot, Khaber, fP, fPi, gP, gPi, Pna1,  Roul,RSU,  SuS, Tna1 `):
 print(``):
else
ezra(args):
fi:
end:

ezra:=proc()
if args=NULL then

print(` ParkingAndTrees.txt: A Maple package for studying Parking Functions and their relations to trees and forests `):
print(`The MAIN procedures are: Cn, Pn, cn, pn, Cna, Pna, cna, pna, Tna, PtoT, TtoP, TnaGeneral, Inv, PairDist, PnzT, PnzF, Pnax, PnPD, SumOfHeight, TotalHeightT, TotalHeightF, ImagePtoT, ImagePtoTall, ImageTtoP, ImageTtoPboth, TnEB, FnEB, CanF, DrawRT, Milon, MilonV, MilonVV, PnaNew, PtoTEB, PtoT0, RPF, RRT, tn, TnaEB, Tn0`):
print(``):

elif nargs=1 and args[1]=Cn then
print(`Cn(n): the set of weakly increasing sequences of length n such that the ith entry of the sequence >=1 and <=i, i.e. weakly increasing parking functions of length n. Try: `):
print(`Cn(5);`):

elif nargs=1 and args[1]=Pn then
print(`Pn(n): the set of parking functions of length n. The elements are obtained by permuting each element in Cn(n).  Try:`):
print(`Pn(4); `):

elif nargs=1 and args[1]=Cna then
print(`Cna(n,a): the set of weakly increasing sequences of length n such that the ith entry of a sequence >=1 and <=a+i-1. Try:`):
print(`Cna(4,3);`):

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

elif nargs=1 and args[1]=cna then
print(`cna(n,a): the number of weakly increasing a-parking functions of length n obtained by recurrence relations.`):
print(`By considering the number of 1’s in the sequence there is a recurrence relation cna(n,a)=add(cna(n-k,a+k-1),k=0..n). Try:`):
print(`cna(4,3);`):

elif nargs=1 and args[1]=pna then
print(`pna(n,a): the number of a-parking functions of length n obtained by recurrence relations.`):
print(`By considering the number of 1’s in the sequence there is a recurrence relation pna(n,a)=add(binomial(n,k)*pna(n-k,a+k-1),k=0..n). Try:`):
print(`pna(4,3);`):

elif nargs=1 and args[1]=Tna then
print(`Tna(n,a): the set of rooted labelled forests on (a+n) vertices with a components, where the roots must be labelled from 1 to a. Try:`):
print(`Tna(4,3);`):

elif nargs=1 and args[1]=PtoT then
print(`PtoT(L,a): map an a-parking function L (as a list) to a forest in Tna(n,a) non-recursively. Try:`):
print(`PtoT([5,8,4,2,1,2,1],2);`):

elif nargs=1 and args[1]=TtoP then
print(`TtoP(T): map a forest T in Tna(n,a) to an a-parking function non-recursively. Try:`):
print(`TtoP({[1,{{1,7},{1,9},{5,9}}],[2,{{2,6},{2,8},{3,4},{3,6}}]});`):

elif nargs=1 and args[1]=TnaGeneral then
print(`TnaGeneral(n,a): the set of forests of rooted labelled trees on (a+n) vertices with a components, where there is no restriction on which vertices can be the roots. Try:`):
print(`Tna(2,2);`):

elif nargs=1 and args[1]=Inv then
print(`Inv(T): input a forest in TnaGeneral(n,a) (if connected then it is a tree, i.e. a=1), output the inversion number of the forest. Try:`):
print(`Inv({[3, {{2, 3}}], [4, {{1, 4}, {1, 5}}]});`):

elif nargs=1 and args[1]=PairDist then
print(`PairDist(T): input a tree, output the total pair distances over all pairs of its vertices. Try:`):
print(`PairDist({[4, {{1, 2}, {1, 3}, {1, 4}}]});`):

elif nargs=1 and args[1]=PnzT then
print(`PnzT(z,n): the sum of z^Inv(T) over all labelled rooted trees on [n]. Try:`):
print(`PnzT(z,4);`):

elif nargs=1 and args[1]=PnzF then
print(`PnzF(z,n): the sum of z^Inv(T) over all labelled rooted forests on [n]. Try:`):
print(`PnzF(z,4);`):

elif nargs=1 and args[1]=Pnax then
print(`Pnax(n,a,x): the sum of x^(sum of all entries in the parking function) over the set of a-parking functions of length n by recurrence relation. Try:`):
print(`Pnax(3,2,x);`):

elif nargs=1 and args[1]=PnPD then
print(`PnPD(z,n): the sum of z^PairDist(T) over all labelled rooted trees on [n]. Try:`):
print(`PnPD(z,5);`):

elif nargs=1 and args[1]=SumOfHeight then
print(`SumOfHeight(T): input a tree or forest, output the sum of heights of all vertices. Try:`):
print(`SumOfHeight({[1, {{1, 4}}], [5, {{2, 5}, {3, 5}}]});`):

elif nargs=1 and args[1]=TotalHeightT then
print(`TotalHeightT(x,n): the sum of x^SumOfHeight(T) over all labelled rooted trees on [n]. Try:`):
print(`TotalHeightT(x,4);`):

elif nargs=1 and args[1]=TotalHeightF then
print(`TotalHeightF(x,n): the sum of x^SumOfHeight(T) over all labelled rooted forests on [n]. Try:`):
print(`TotalHeightF(x,4);`):

elif nargs=1 and args[1]=ImagePtoT then
print(`ImagePtoT(n,a): output the image of all a-parking functions of length n under the map PtoT. Try:`):
print(`ImagePtoT(3,2);`):

elif nargs=1 and args[1]=ImagePtoTall then
print(`ImagePtoTall(n,a): output the image of all a-parking functions of length n under the map PtoT and the pre-image together. Try:`):
print(`ImagePtoTall(3,2);`):


elif nargs=1 and args[1]=ImagePtoTboth then
print(`ImagePtoTboth(n): output the image of all 1-parking functions of length n under the map PtoT and PtoRT0(P),  Try:`):
print(`ImagePtoTboth(3);`):

elif nargs=1 and args[1]=ImageTtoP then
print(`ImageTtoP(n,a): output the image of all labelled rooted forests on [a+n] with a components where [a] are the a roots under the map TtoP. Try:`):
print(`ImageTtoP(3,2);`):

elif nargs=1 and args[1]=ImageTtoPboth then
print(`ImageTtoPboth(n,a): output the image of all labelled rooted forests on [a+n] with a components where [a] are the a roots under the map TtoP and the pre-image together. Try:`):
print(`ImageTtoPboth(3,2);`):

elif nargs=1 and args[1]=TnEB then
print(`TnEB(n): the set of labelled rooted trees on [n] by considering which vertex is the root and using Tna1(Emps, Bos). Try:`):
print(`TnEB(4);`):

elif nargs=1 and args[1]=FnEB then
print(` FnEB(n): the set of labelled rooted forests on [n] by considering which vertices are the roots and using Tna1(Emps, Bos). Try:`):
print(`FnEB(4);`):

elif nops([args])=1 and op(1,[args])=CanF then
print(`CanF(T,r): the canonical form of a rooted tree T with root r. Try: `):
print(`CanF({{0,1},{0,2},{0,3},{1,4}},0);`):

elif nops([args])=1 and op(1,[args])=CheckfP then
print(`CheckfP(n,a): checks that fPi is the inverse of fP (and vice versa), for all a-parking functions of length n. Try: `):
print(`CheckfP(4,1); `):

elif nops([args])=1 and op(1,[args])=CheckgP then
print(`CheckgP(Emps,Bos): checks that gPi is the inverse of gP (and vice versa), for all forests of rooted trees with set of roots Bos and`):
print(`set of NonRoots Emps. Try: `):
print(`CheckgP({1,2},{3,4}); `):

elif nops([args])=1 and op(1,[args])=CheckPtoTEB then
print(`CheckPtoTEB(n,a): Checks PtoTEB for a set of Employees with n elements and a set of bosses with a elements. Try:`):
print(`CheckPtoTEB(4,1);`):

elif nops([args])=1 and op(1,[args])=Cna1 then
print(`Cna1(n,a): the set of sorted a-parking functions. Try:`):
print(`Cna1(3,1);`):

elif nops([args])=1 and op(1,[args])=Conv1 then
print(`Conv1(F): converts a forest or rooted trees F into the format [SetOfNonRoots, SetOfRoots,SetOfEdges]. Try:`):
print(`Conv1({[1,{{1,2},{1,3}}]);`):

elif nops([args])=1 and op(1,[args])=Dorot then
print(`Dorot(T): inputs a tree, in canonical form, and outputs, in order, the successive generations. Try:`):
print(`Dorot(CanF({{0,1},{0,2},{0,3},{1,4}},0));`):

elif nops([args])=1 and op(1,[args])=DrawRT then
print(`DrawRT(T,r): draws the rooted tree T with root r. Try: `):
print(` DrawRT({{0,1},{0,2},{0,3}},0); `):

elif nops([args])=1 and op(1,[args])=fP then
print(`fP(P): inputs a pair  PP=[a,P] where a is >=1 and P is a-parking function of length n=nops(P). Outputs the pair [S,[a+k-1,P1]] where`):
print(`S is the subset of {1, ...n} of i's such that P[i]=1 and P1 is the a+nops(S)-1 parking function P1 of length n-nops(S) obtained by`):
print(`subtracting one from each entry. Try:`):
print(`fP([1,[1,2,3,4]]); `):

elif nops([args])=1 and op(1,[args])=fPi then
print(`fPi(sPP): inputs a pair sPP=[S,[a1,P1]] where S is a set of integers, subset of [1,nops(P1)+nops(S)] and outputs`):
print(`a pair [a,P] where P is an a-parking function. It is the inverse of fP(P). Try:`):
print(`fPi(fP([1,[1,2,3,4]])) ;`):

elif nops([args])=1 and op(1,[args])=gP then
print(`gP(F): inputs a forest F (in new format [Emps,Bos,SetOfEdges]) and outputs the pair [S,F1] where`):
print(`S is the subset of Emps connected to the smallest member of Bos,let's call it b1, and F1 is the tree `):
print(`[Emps minus S,Bos minus {b1} union S,SetOfEdges without b1]. Try:`):
print(`gP([{2,3,4},{1},{{1,2},{1,3},{1,4}}]);`):

elif nops([args])=1 and op(1,[args])=gPi then
print(`gPi(Emps,Bos,P): inputs a set of Employees,Emps, a set of Bosses Bos, and a forest of rooted trees in new format, outputs the inverse of gP. Try: `):
print(` gPi({2,3,4},{1}, gP([{2,3,4},{1},{{1,2},{1,3},{1,4}}]));  `):

elif nops([args])=1 and op(1,[args])=Khaber then
print(`Khaber(S,L): combines a set S and a list L by creating a new list of size nops(S)+nops(L) such that`):
print(`the positions that belong to S are all 1 and the other positions are those of L +1, in that order. Try: `):
print(`Khaber({1,3},[1,2]);`):

elif nops([args])=1 and op(1,[args])=Milon then
print(`Milon(n): all the images of the 1-parking functions of length n as rooted trees with 0 the root. Try:`):
print(`Milon(4);`):

elif nops([args])=1 and op(1,[args])=MilonV then
print(`MilonV(n): all the images of the 1-parking functions of length n as rooted trees with 0 the root, arranged by sorting, verbose version. Try:`):
print(`MilonV(3);`):


elif nops([args])=1 and op(1,[args])=MilonVV then
print(`MilonVV(n): all the images of the 1-parking functions of length n as rooted trees with 0 the root, arranged by sorting, verbose version. `):
print(` Also draws them. Try:`):
print(`MilonVV(3);`):

elif nops([args])=1 and op(1,[args])=PnaNew then
print(`PnaNew(n,a): the set of a-parking functions in the format [a, ParkingFunction]. Try:`):
print(`PnaNew(3,1);`):

elif nops([args])=1 and op(1,[args])=Pna1 then
print(`Pna1(n,a): the set of a-parking functions. Try:`):
print(`Pna1(3,1);`):

elif nops([args])=1 and op(1,[args])=PtoTEB then
print(`PtoTEB(P,Emps,Bos): inputs an a:=nops(Bos)-parking function and outputs the corresponding forest of rooted treese with a components.`):
print(`with a prescribed set of Employers and Bosses, Emps and Bos resepctively`):
print(`Try: `):
print(`PtoTEB([2,[2,2,4]],{3,4,5},{1,2});`):

elif nops([args])=1 and op(1,[args])=PtoT0 then
print(`PtoT0(P): given a 1-parking function P, outputs PtoT([1,P],{0},{1, ..., n}), where n=nops(P)). Try:`):
print(`PtoT0([4,3,1,1]); `);


elif nops([args])=1 and op(1,[args])=Roul then
print(`Roul(L): inputs a list L and outputs i from 1 to nops(L) with prob. L[i]/convert(L,`+`); Try: `):
print( `Roul([1,2,1]); `):

elif nops([args])=1 and op(1,[args])=RPF then
print(`RPF(n,a): a random a-parking function of length n. Try:`):
print(`RPF(5,1);`):

elif nops([args])=1 and op(1,[args])=RRT then
print(`RRT(n): a random labelled rooted tree with root 0 and non-roots labelled {1, ..., n}. Try:`):
print(`RRT(10);`):

elif nops([args])=1 and op(1,[args])=RSU then
print(`RSU(n,k): a random k-subset  of {1, ..., n}. Try: RSU(5,3);`):

elif nops([args])=1 and op(1,[args])=SuS then
print(`SuS(S0,U): inputs a subset S0 of {1, ..., nops(U)} and outputs the corresponding subset of U. Try`):
print(`SuS({1,3},{4,5,7,8});`):

elif nops([args])=1 and op(1,[args])=tn then
print(`tn(n,a) the number of forests of rooted trees  where each component has the form [root, SetOfEdges]`):
print(`n employees and a bosses, i.e. the number of ordered labelled forests`):

elif nops([args])=1 and op(1,[args])=Tna1 then
print(`Tna1(Emps,Bos) The set of sets of rooted trees where the set of  roots is Bos and the set of non-root vertices is Emps. Try`):
print(`Tna1({1,2},{3,4,5});`):

elif nops([args])=1 and op(1,[args])=TnaEB then
print(`TnaEB(Emps,Bos) The set of sets of rooted trees where the set of  roots is Bos and the set of non-root vertices is Emps. New format. Try`):
print(`TnaEB({1,2},{3,4,5});`):

elif nops([args])=1 and op(1,[args])=Tn0 then
print(`Tn0(n) The set of rooted trees on {0,...,n} rooted at 0. Try:`):
print(`Tn0(3);`):

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


#Cn(n): the set of weakly increasing sequences of length n such that the ith entry of the sequence >=1 and <=i.
Cn:=proc(n) local S,L,l,j:
option remember:
if n=0 then
  return {[]}:
fi:
if n=1 then
  return {[1]}:
fi:
S:={}:
L:=Cn(n-1):
for l in L do
  S:=S union {seq([op(l),j],j=l[-1]..n)}:
od:
S:
end:

#numCn(n): number of elements in Cn(n) by using nops.
numCn:=proc(n)
nops(Cn(n)):
end:


#Pn(n): the set of length n sequences that are obtained by permuting each element in Cn(n).
Pn:=proc(n) local S,L,l:
S:={}:
L:=Cn(n):
for l in L do 
  S:=S union {op(permute(l))}:
od:
end:

#numPn(n): number of elements in Pn(n),i.e. parking function of length n by using nops.
numPn:=proc(n)
nops(Pn(n)):
end:


#cnm(n,m): number of weakly increasing sequences of length n such that the ith entry <= i and the last entry is m.
cnm:=proc(n,m) local k:
option remember:
if m>n then 
  return 0:
fi:
if m<=0 then 
  return 0:
fi:
if m=1 and n>=1 then
  return 1:
fi:
add(cnm(n-1,k),k=1..m): 
end:

#cn(n): number of weakly increasing sequences of length n such that the ith entry <= i without using nops().
cn:=proc(n) local m:
add(cnm(n,m),m=1..n):
end:


#pnm(n,m): number of sequences of length n such that for any i in [1,n] the number of elements >i is at most n-i and the largest element is m.
pnm:=proc(n,m) local k,i:
option remember:
if m>n then
  return 0:
fi:
if m<=0 then
  return 0:
fi:
if m=1 and n>=1 then
  return 1:
fi:
add(binomial(n,k)*add(pnm(n-k,i),i=1..m-1),k=1..n+1-m): #depends on the number of m in the sequence
end:

#pn(n): number of parking functions of length n.
pn:=proc(n) local m:
add(pnm(n,m),m=1..n):
end:


#Cna(n,a): the set of weakly increasing sequences of length n such that the ith entry of a sequence >=1 and <=a+i-1.
Cna:=proc(n,a) local i,S,L,l,j:
option remember:
if n<0 or a<1 then
  return FAIL:
fi:
if n=0 then
  return {[]}:
fi:
if n=1 then 
  return {seq([i],i=1..a)}:
fi:
S:={}:
L:=Cna(n-1,a):
for l in L do
  S:=S union {seq([op(l),j],j=l[-1]..n+a-1)}:
od:
S:
end:
  
#cna(n,a): number of weakly increasing sequences of length n such that the ith entry of a sequence >=1 and <=a+i-1. By considering the number of 1’s in the sequence there is a recurrence relation cna(n,a)=add(cna(n-k,a+k-1),k=0..n).
cna:=proc(n,a) local k:
option remember:
if n=0 then
  return 1:
fi:
if n>=1 and a<=0 then
  return 0:
fi:
add(cna(n-k,a+k-1),k=0..n):
end:


#Pna(n,a): the set of sequences of length n that are obtained by permuting each element of Cna(n,a).
Pna:=proc(n,a) local S,L,l:
S:={}:
L:=Cna(n,a):
for l in L do 
  S:=S union {op(permute(l))}:
od:
end:

#pna(n,a): number of sequences of length n that are obtained by permuting each element of Cna(n,a).
pna:=proc(n,a) local k:
option remember:
if n=0 then
  return 1:
fi:
if n>=1 and a<=0 then
  return 0:
fi:
add(binomial(n,k)*pna(n-k,a+k-1),k=0..n):
end:
#pna(n,a) = a*(a+n)^(n-1)


#TnaG(n,a): a set of list where the first element in the list is the set of forests of rooted labelled trees on (a+n) vertices with a components, where the roots must be labelled from 1 to a. The second element in the list is a list of cardinality a whose ith element is a set of numbers such that the vertex labelled by this number is in the tree with root i (1<=i<=a).
TnaG:=proc(n,a) local i,S,T,t,F,L,l,j,k,newF,newL:
if a<1 or n<0 then 
  return FAIL:
fi:
if n=0 then
  return {[{seq([i,{}],i=1..a)},[seq({i},i=1..a)]]}:
fi:

S:={}:
T:=TnaG(n-1,a):

for j from a+1 to a+n do
  for t in T do
  t:=subs(j=a+n,t):
  F:=t[1]:
  L:=t[2]:
    for k from 1 to a do
      if k<>1 and k<>a then
        newL:=[op(1..k-1,L),L[k] union {j},op(k+1..a,L)]:
        for l in L[k] do
          newF:={op(1..k-1,F),[F[k][1],F[k][2] union {{l,j}}],op(k+1..a,F)}:
          S:=S union {[newF,newL]}:
        od:
      elif k=1 then
        newL:=[L[k] union {j},op(k+1..a,L)]:
        for l in L[k] do
          newF:={[F[k][1],F[k][2] union {{l,j}}],op(k+1..a,F)}:
          S:=S union {[newF,newL]}:
        od:
      elif k=a then
        newL:=[op(1..k-1,L),L[k] union {j}]:
        for l in L[k] do
          newF:={op(1..k-1,F),[F[k][1],F[k][2] union {{l,j}}]}:
          S:=S union {[newF,newL]}:
        od:
      fi:
    od:
  od:
od:
S:
end:
        
#Tna(n,a): the set of forests of rooted labelled trees on (a+n) vertices with a components, where the roots must be labelled from 1 to a.
Tna:=proc(n,a) local F:
{seq(F[1],F in TnaG(n,a))}:
end:


#PtoS(L): map an a-parking function to a list of set by considering the index of 1 in the function recursively. L must be an a-parking function.
PtoS1step:=proc(L) local n,i,S,L1,nonOneIndices:
if L=[] then
  return [{},[]]:
fi:
n:=nops(L):
S:={seq(`if`(L[i]<>1,NULL,i),i=1..n)}:
nonOneIndices:=[seq(`if`(L[i]=1,NULL,i),i=1..n)]:
L1:=L[nonOneIndices]:
L1:=L1+[-1$nops(L1)]:
[S,L1]:
end:

PtoS:=proc(L) local T,LS:
T:=L:
LS:=[]:
while PtoS1step(T)<>[{},[]] do
  LS:=[op(LS),PtoS1step(T)[1]]:
  T:=PtoS1step(T)[2]:
od:
LS:
end:


#StoP(S): map a list of set to an a-parking function, which is the inverse of the above PtoS(L) procedure.
StoP:=proc(S) local n,L,i,j:
n:=nops(S):
L:=[]:
for i from n to 1 by -1 do
  if S[i]={} then
    L:=L+[1$nops(L)]:
  else
    L:=L+[1$nops(L)]:
    for j in S[i] do
      if j=1 then 
        L:=[1,op(L)]:
      elif j=1+nops(L) then
        L:=[op(L),1]:
      else
        L:=[op(1..j-1,L),1,op(j..nops(L),L)]:
      fi:
    od:
  fi:
od:
L:
end:
     

#TtoS(T): map a forest in Tna(n,a) to a list of sets by considering the neighbors of vertex 1 recursively. #not sure about the index of the new forest after each step, might be wrong

#Given a vertex v and a set of edges S, find those edges which can be in a tree with v.
ConsTree:=proc(v,S) local S1,C,edge,child,T:
if not v in `union`(seq(edge, edge in S)) then
  return {}:
fi:
C:={}:
T:={}:
S1:=S:
for edge in S1 do
  if v in edge then
    T:=T union {edge}:
    C:=C union edge minus {v}:
    S1:=S1 minus {edge}:
  fi:
od:
return (`union`(seq(ConsTree(child,S1),child in C)) union T):
end:

TtoS1step:=proc(T) local Eset,edge,S,newF,s,newT,candidate,V,v,tree,A,B,i:
if T={} then
  return 0:
fi:
Eset:=T[1][2]:
S:={}:
for edge in Eset do
  if edge[1]=1 then 
    S:=S union {edge[2]}:
  fi:
od:
candidate:={}: #all edges in tree 1 without vertex 1
for edge in Eset do
  if not 1 in edge then
    candidate:=candidate union {edge}:
  fi:
od:
newF:={op(2..nops(T),T)}:
for s in S do
  newT:=[s,ConsTree(s,candidate)]:
  newF:=newF union {newT}:
od:
V:=`union`(seq(tree[2],tree in T)):
V:=`union`(seq(v,v in V)) union {seq(i,i=1..nops(T))}: #the set of all vertices
A:=[seq(i,i=1..nops(V)-1)]:
B:=[seq(i,i=2..nops(T)),op(S),op(V minus S minus {seq(i,i=1..nops(T))})]:
newF:=subs([seq(B[i]=A[i],i=1..nops(A))],newF):
[S,newF]:
end:

TtoS:=proc(T) local S,newT:
newT:=T:
S:=[]:
while TtoS1step(newT)<>0 do
  S:=[op(S),TtoS1step(newT)[1]]:
  newT:=TtoS1step(newT)[2]:
od:
S:
end:


#StoT(S): map a list of sets to a forest in Tna(n,a).    #unfinished 
#StoT:=proc(S) local F,n,i,numT:
#n:=nops(S):
#F:={}:
#for i from n to 1 by -1 do
  #numT:=nops(F):
  



#PtoT(L,a): map an a-parking function to a forest in Tna(n,a).
PtoT:=proc(L,a) local F,i,n,wiL,A,B,V,v,j,L1,V1,V2,T:
F:={seq([i,{}],i=1..a)}: #F is the forest to output
n:=nops(L): 
V:=[seq(i,i=a+1..a+n)]: #vertices which are not roots
wiL:=sort(L):
A:=table([seq(V[i]=wiL[i],i=1..n)]):
B:=table([seq(V[i]=L[i],i=1..n)]):
L1:=[seq({i},i=1..a)]:
for v in V do
  for j from 1 to a do
    if A[v] in L1[j] then
      L1[j]:=L1[j] union {v}:
      T:=F[j]:
      T[2]:=T[2] union {{v,A[v]}}:
      F:=F minus {F[j]} union {T}:
      break:
    fi:
  od:
od:
V1:=V:
V2:=[]:
for i from 1 to n do
  for j from 1 to nops(V1) do
    if B[V1[j]]=wiL[i] then
      V2:=[op(V2),V1[j]]:
      if j=1 then
        V1:=[op(2..nops(V1),V1)]:
      elif j=nops(V1) then
        V1:=[op(1..nops(V1)-1,V1)]:
      else
        V1:=[op(1..j-1,V1),op(j+1..nops(V1),V1)]:
      fi:
      break:
    fi:
  od:
od:
F:=subs([seq(V[i]=V2[i],i=1..n)],F):
end:


#TtoP(T): map a forest in Tna(n,a) to an a-parking function.
Depth:=proc(T,n,a) #input a forest in Tna(n,a) output a list L of length a+n s.t. that L[i] is the depth of i in the forest. Set Depth(1)=0.
local L,large,Edges,edge,i,Known,New:
large:=2*(a+n)+100:
L:=[0$a,large$n]:
Edges:=`union`(seq(T[i][2],i=1..a)):
Known:={seq(i,i=1..a)}:
for i from 1 to n do
  if large in L then
    New:={}:
    for edge in Edges do
      if edge intersect Known <> {} then
        if edge[1] in Known and L[edge[2]]=large then
          L[edge[2]]:=L[edge[1]]+1:
          New:=New union {edge[2]}:
        elif edge[2] in Known and L[edge[1]]=large then
          L[edge[1]]:=L[edge[2]]+1:
          New:=New union {edge[1]}:
        fi:
      fi:
    od:
    Known:=New:
  fi:
od:
L:
end:

Parent:=proc(T,n,a) #input a forest in Tna(n,a) output a list L of length a+n s.t. that L[i] is the parent of i in the forest. Set Parent(1..a)=0.
local L,i,j,Edges,D,d,pd,Candidate,cand:
L:=[0$(a+n)]:
Edges:=`union`(seq(T[i][2],i=1..a)):
D:=Depth(T,n,a):
for j from a+1 to a+n do
  d:=D[j]:
  pd:=d-1:
  Candidate:={seq(`if`(D[i]<>pd,NULL,i),i=1 .. a+n)}:
  for cand in Candidate do
    if {j,cand} in Edges then
      L[j]:=cand:
      break:
    fi:
  od:
od:
L:
end:

Index:=proc(T,n,a) #input a forest in Tna(n,a) output a list L of length a+n s.t. that L[i] is the vertex at the same position in the basic forest. Set index(1..a)=1..a.
local L,i,j,m,k,l,D,P,S,Consider,c:
L:=[seq(i,i=1..a),0$n]:
D:=Depth(T,n,a):
P:=Parent(T,n,a):
m:=max(D):
S:={seq(i,i=a+1..a+n)}:
for j from 1 to m do
  Consider:=[seq(`if`(D[i]<>j,NULL,i),i=1 .. a+n)]:
  if nops(Consider)=0 then
    break:
  fi:
  if nops(Consider)=1 then
    L[Consider[1]]:=S[1]:
    S:=S minus {S[1]}:
    next:
  fi:
  for k from 1 to nops(Consider)-1 do
    for l from nops(Consider) to k+1 by -1 do
      if L[P[Consider[l]]]<L[P[Consider[l-1]]] or (L[P[Consider[l]]]=L[P[Consider[l-1]]] and Consider[l]<Consider[l-1]) then
        Consider[l],Consider[l-1]:=Consider[l-1],Consider[l]:
      fi:
    od:
  od:
  for c in Consider do
    L[c]:=S[1]:
    S:=S minus {S[1]}:
  od:
od:
L:
end:

TtoP:=proc(T) local a,n,i,V,L,P,Ind:
a:=nops(T):
V:=`union`(seq(T[i][2],i=1..a)):
V:=`union`(seq(i,i in V)) union {seq(i,i=1..a)}:
n:=nops(V)-a:
L:=[seq(i,i=1..a),0$n]:
P:=Parent(T,n,a):
Ind:=Index(T,n,a):
for i from a+1 to a+n do
  L[i]:=Ind[P[i]]:
od:
L:=L[a+1..a+n]:
end:


#TnaGeneral(n,a): the set of forests of rooted labelled trees on (a+n) vertices with a components, where there is no restriction on which vertices can be the roots.
TnaGeneral:=proc(n,a) local L,S,F,l,f,A,i:
F:=Tna(n,a):
L:=permute(n+a):
S:={}:
A:=[seq(i,i=1..n+a)]:
for f in F do
  for l in L do
    S:=S union {subs({seq(A[i]=l[i],i=1..n+a)},f)}:
  od:
od:
S:
end:

#CheckTtoP(n,a): check the image of all t’s in Tna(n,a) under TtoP is indeed Pna(n,a). 
CheckTtoP:=proc(n,a) local t:
evalb({seq(TtoP(t), t in Tna(n, a))} = Pna(n, a)):
end:

#CheckPtoT(n,a): check the image of all p’s in Pna(n,a) under PtoT is indeed Tna(n,a). 
CheckPtoT:=proc(n,a) local p:
evalb({seq(PtoT(p, a), p in Pna(n, a))} = Tna(n, a)):
end:

#CheckPinverse(n,a): check whether TtoP(PtoT(p))=p for all p’s in Pna(n,a).
CheckPinverse:=proc(n,a) local p:
op({seq(evalb(TtoP(PtoT(p, a)) = p), p in Pna(n, a))}):
end:

#CheckTinverse(n,a): check whether PtoT(TtoP(t))=t for all t’s in Tna(n,a).
CheckTinverse:=proc(n,a) local t:
op({seq(evalb(PtoT(TtoP(t), a) = t), t in Tna(n, a))}):
end:


#Inv(T): input a forest in TnaGeneral(n,a) (if connected then it is a tree, i.e. a=1), output the inversion number of the forest.

depth:=proc(T,n,a) #input a forest in TnaGeneral(n,a) output a list L of length a+n s.t. that L[i] is the depth of i in the forest. Set depth(roots)=0.
local L,large,Edges,edge,i,Known,New:
large:=2*(a+n)+100:
L:=[large$(n+a)]:
for i from 1 to a do
  L[T[i][1]]:=0:
od:
Edges:=`union`(seq(T[i][2],i=1..a)):
Known:={seq(T[i][1],i=1..a)}:
for i from 1 to n do
  if large in L then
    New:={}:
    for edge in Edges do
      if edge intersect Known <> {} then
        if edge[1] in Known and L[edge[2]]=large then
          L[edge[2]]:=L[edge[1]]+1:
          New:=New union {edge[2]}:
        elif edge[2] in Known and L[edge[1]]=large then
          L[edge[1]]:=L[edge[2]]+1:
          New:=New union {edge[1]}:
        fi:
      fi:
    od:
    Known:=New:
  fi:
od:
L:
end:

parent:=proc(T,n,a) #input a forest in TnaGeneral(n,a) output a list L of length a+n s.t. that L[i] is the parent of i in the forest. Set parent(roots)=0.
local L,i,j,J,Edges,D,d,pd,Candidate,cand:
L:=[0$(a+n)]:
Edges:=`union`(seq(T[i][2],i=1..a)):
D:=depth(T,n,a):
J:={seq(i,i=1..a+n)} minus {seq(T[i][1],i=1..a)}:
for j in J do
  d:=D[j]:
  pd:=d-1:
  Candidate:={seq(`if`(D[i]<>pd,NULL,i),i=1 .. a+n)}:
  for cand in Candidate do
    if {j,cand} in Edges then
      L[j]:=cand:
      break:
    fi:
  od:
od:
L:
end:

ancestor:=proc(T,n,a) #input a forest in TnaGeneral(n,a) output a list L of length a+n s.t. that L[i] is the set of ancestors of i in the forest. Set ancestor(roots)={}.
local L,D,P,m,i,j,Temp,temp:
L:=[{}$(a+n)]:
D:=depth(T,n,a):
P:=parent(T,n,a):
m:=max(D):
for j from 1 to m do
  Temp:={seq(`if`(D[i]<>j,NULL,i),i=1 .. a+n)}:
  for temp in Temp do
    L[temp]:=L[P[temp]] union {P[temp]}:
  od:
od:
L:
end:

Inv:=proc(T) local a,n,V,i,co,j,A:
a:=nops(T):
V:=`union`(seq(T[i][2],i=1..a)):
V:=`union`(seq(i,i in V)) union {seq(T[i][1],i=1..a)}:
n:=nops(V)-a:
co:=0:
A:=ancestor(T,n,a):
for i in V do
  for j in A[i] do
    if j>i then
      co:=co+1:
    fi:
  od:
od:
co:
end:


#PairDist(T): input a tree, output the total pair distances over all pair of its vertices. 

pairdist:=proc(T,x,y) #Calculate the pair distance of (x,y) on the tree T.
local a,n,V,i,A,D,Common,M,c,CA:
a:=nops(T):
V:=`union`(seq(T[i][2],i=1..a)):
V:=`union`(seq(i,i in V)) union {seq(T[i][1],i=1..a)}:
n:=nops(V)-a:
if a<>1 or not x in V or not y in V then
  return fail:
fi:
if x=y then
  return 0:
fi:
A:=ancestor(T,n,a):
D:=depth(T,n,a):
if x in A[y] or y in A[x] then
  return abs(D[x]-D[y]):
fi:
Common:=A[x] intersect A[y]:
M:=max({seq(D[c],c in Common)}):
CA:=op({seq(`if`(D[c]<>M,NULL,c),c in Common)}):
return D[x]+D[y]-2*D[CA]:
end:

PairDist:=proc(T) local i,j,a,n,V,L,s:
a:=nops(T):
V:=`union`(seq(T[i][2],i=1..a)):
V:=`union`(seq(i,i in V)) union {seq(T[i][1],i=1..a)}:
n:=nops(V)-a:
s:=0:
if a<>1 then 
  return fail:
fi:
L:=[op(V)]:
for i from 1 to nops(L)-1 do
  for j from i+1 to nops(L) do
    s:=s+pairdist(T,L[i],L[j]):
  od:
od:
s:
end:

#PnzT(z,n): the sum of z^Inv(T) over all labelled rooted trees on [n].
PnzT:=proc(z,n) local t:
add(z^Inv(t), t in TnEB(n)):
end:

#PnzT1(z,n): the sum of z^Inv(T) over all labelled rooted trees on [n].
PnzT1:=proc(z,n) local t:
add(z^Inv(t), t in TnaGeneral(n-1,1)):
end:


#PnzF(z,n): the sum of z^Inv(F) over all labelled rooted forests on [n].
PnzF:=proc(z,n) local f:
add(z^Inv(f), f in FnEB(n)): 
end:

#PnzF1(z,n): the sum of z^Inv(F) over all labelled rooted forests on [n].
PnzF1:=proc(z,n) local f:
add(z^Inv(f), f in Tna(n,1)): #add a root 1 to the forest to be a tree in Tna(n,1).
end:


#Pnax(n,a,x): the sum of x^(sum of all entries in the parking function) over the set of a-parking functions of length n by recurrence relation.
Pnax:=proc(n,a,x) local k:
option remember:
if n=0 then
  return 1:
fi:
if n>0 and a=0 then
  return 0:
fi:
return expand(x^n*add(binomial(n,k)*Pnax(n-k,a+k-1,x),k=0..n)):
end:

#PnPD(z,n): the sum of z^PairDist(T) over all labelled rooted trees on [n].
PnPD:=proc(z,n) local t:
add(z^PairDist(t), t in TnEB(n)):
end:

#SumOfHeight(T): input a tree or forest, output the sum of heights of all vertices.
SumOfHeight:=proc(T) local a,V,n,D,i:
a:=nops(T):
V:=`union`(seq(T[i][2],i=1..a)):
V:=`union`(seq(i,i in V)) union {seq(T[i][1],i=1..a)}:
n:=nops(V)-a:
D:=depth(T,n,a):
convert(D,`+`):
end:


#TotalHeightT(x,n): the sum of x^SumOfHeight(T) over all labelled rooted trees on [n].
TotalHeightT:=proc(x,n) local T:
add(x^SumOfHeight(T), T in TnaGeneral(n-1,1)):
end:

#TotalHeightF(x,n): the sum of x^SumOfHeight(T) over all labelled rooted forests on [n].
TotalHeightF:=proc(x,n) local F:
add(x^(SumOfHeight(F)-n), F in Tna(n,1)):
end:


#ImagePtoT(n,a): output the image of all a-parking functions of length n under the map PtoT.
ImagePtoT:=proc(n,a) local P:
for P in Pna(n,a) do
  print(PtoT(P,a)):
od:
end:

#ImagePtoTall(n,a): output the image of all a-parking functions of length n under the map PtoT and the pre-image together.
ImagePtoTall:=proc(n,a) local P:
for P in Pna(n,a) do
  print(P, PtoT(P,a)):
od:
end:


#ImagePtoTboth(n): output the image of all a-parking functions of length n under the map PtoT and the pre-image together.
ImagePtoTboth:=proc(n) local P,i:
for P in Pna(n,1) do
  print(P, subs({seq(i=i-1,i=1..n+1)},PtoT(P,1)[1][2]),PtoT0(P)):
od:
end:

#ImageTtoP(n,a): output the image of all labelled rooted forests on [a+n] with a components where [a] are the a roots under the map TtoP.
ImageTtoP:=proc(n,a) local T:
for T in Tna(n,a) do
  print(TtoP(T)):
od:
end:

#ImageTtoPboth(n,a): output the image of all labelled rooted forests on [a+n] with a components where [a] are the a roots under the map TtoP and the pre-image together.
ImageTtoPboth:=proc(n,a) local T:
for T in Tna(n,a) do
  print(T, TtoP(T)):
od:
end:


#TnEB(n): the set of labelled rooted trees on [n] by considering which vertex is the root and using Tna1(Emps, Bos). This procedure should be faster than TnaGeneral(n,a). 
TnEB:=proc(n) local S,i,N:
S:={}:
N:={seq(i,i=1..n)}:
for i from 1 to n do
  S:=S union Tna1(N minus {i}, {i}):
od:
S:
end:


#FnEB(n): the set of labelled rooted forests on [n] by considering which vertices are the roots and using Tna1(Emps, Bos). This procedure should be faster than TnaGeneral(n,a). 
FnEB:=proc(n) local N,i,PO,P,S:
N:={seq(i,i=1..n)}:
PO:=powerset(N) minus {{}}:
S:={}:
for P in PO do
  S:=S union Tna1(N minus P, P):
od:
S:
end:


###### The procedures below are from ETZIM.txt written by Doron Zeilberger ######


#tn(n,a) the number of forests of rooted trees  where each component has the form [root, SetOfEdges]
#n employees and a bosses, i.e. the number of ordered labelled forests
tn:=proc(n,a) local k:
option remember:


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

if a<=0 then
 RETURN(0):
fi:

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


#Tna1(Emps,Bos) The set of sets of rooted trees where the set of  roots is Bos and the set of non-root vertices is Emps. Try
#Tna1({1,2},{3,4,5});
Tna1:=proc(Emps,Bos) local k,b,gu,mu,lu,lu1,mu1,mu1a,mu11, FiT:
option remember:

if Bos={} then
 RETURN({}):
fi:

if Emps={} then
 RETURN({{seq([b,{}], b in Bos)}}):
fi:

b:=min(op(Bos)):

mu:=Tna1(Emps,Bos minus {b}):
gu:={seq({[b,{}]} union mu1, mu1 in mu)}:


for k from 1 to nops(Emps) do

 lu:=choose(Emps,k):

 for lu1 in lu do


    mu:=Tna1(Emps minus lu1, Bos minus {b} union lu1):


    for mu1 in mu do

   
       mu1a:=mu1:
       FiT:={}:
       for mu11 in mu1 do


        if member(mu11[1],lu1) then
         mu1a:=mu1a minus {mu11}:

         FiT:=FiT union (mu11[2] union {{b,mu11[1]}}):

          fi:
           od:
 gu:=gu union {mu1a union {[b,FiT]}}:    
od:

od:
od:
gu:
end:

          
#Conv1(F): converts a forest or rooted trees F into the format [SetOfNonRoots, SetOfRoots,SetOfEdges]. Try:
#Conv1({[1,{{1,2},{1,3}}]);
Conv1:=proc(F) local Bo,Edg,edg,ver,i:
 Bo:={seq(F[i][1],i=1..nops(F))}:
 Edg:={seq(op(F[i][2]),i=1..nops(F))}:
 ver:={seq(op(edg),edg in Edg)}:
 [ver minus Bo,Bo,Edg]:
end:

#TnaEB(Emps,Bos) The set of sets of rooted trees where the set of  roots is Bos and the set of non-root vertices is Emps. Try
#TnaEB({1,2},{3,4,5});
TnaEB:=proc(Emps,Bos) local gu,gu1:
option remember:
gu:=Tna1(Emps,Bos):
{seq(Conv1(gu1), gu1 in gu)}:
end:
      

 

#Tn0(n) The set of rooted trees on {0,...,n} rooted at 0. Try
#Tn0(3);
Tn0:=proc(n) local i,gu,gu1:
option remember:
gu:=Tna1({seq(i,i=1..n)},{0}):

{seq(gu1[1][2],gu1 in gu)}:
end:

#Khaber(S,L): combines a set S and a list L by creating a new list of size nops(S)+nops(L) such that
#the positions that belong to S are all 1 and the other positions are those of L +1, in that order. Try
#Khaber({1,3},[1,2]);
Khaber:=proc(S,L) local gu,n,L1,i:
n:=nops(S)+nops(L):
L1:=L:

gu:=[]:

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

gu:
end:


#Pna1(n,a): the set of a-parking functions. Try:
#Pna1(3,1);
Pna1:=proc(n,a) local gu,k,lu1,lu,mu,mu1:
option remember:

if a<1 then
 RETURN({}):
fi:

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

gu:={}:

for k from 0 to n do
lu:=choose(n,k):
 mu:=Pna1(n-k,a+k-1):

for lu1 in lu do
 for mu1 in mu do
  gu:=gu union {Khaber(lu1,mu1)}:
 od:
od:
od:
gu:

end:

#Cna1(n,a): the set of sorted a-parking functions. Try:
#Cna1(3,1);
Cna1:=proc(n,a) local gu,k,mu,mu1,i1:
option remember:

if a<1 then
 RETURN({}):
fi:

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

gu:={}:

for k from 0 to n do
 mu:=Cna1(n-k,a+k-1):
  gu:=gu union {seq([1$k,seq(mu1[i1]+1,i1=1..n-k)],mu1 in mu)}:
od:

gu:

end:



#PnaNew(n,a): the set of a-parking functions. New format Try:
#PnaNew(3,1);
PnaNew:=proc(n,a) local gu,gu1:
option remember:
gu:=Pna1(n,a):
{seq([a,gu1],gu1 in gu)}:
end:


#fP(PP): inputs a pair  PP=[a,P] where a is >=1 and P is a-parking function of length n=nops(P). Outputs the pair [S,[a+k-1,P1]] where
#S is the subset of {1, ...n} of i's such that P[i]=1 and P1 is the a+nops(S)-1 parking function P1 of length n-nops(S) obtained by
#subtracting one from each entry. Try
#fP([1,[1,2,3,4]]);
fP:=proc(PP) local S,P,P1,i,a,n,Ps:
a:=PP[1]:
P:=PP[2]:
n:=nops(P):

Ps:=sort(P):
if {seq(evalb(Ps[i]<=a+i-1),i=1..n)}<>{true} then
  RETURN(FAIL):
fi:

S:=[]:
P1:=[]:

for i from 1 to n do
 if P[i]=1 then
  S:=[op(S),i]:
 else
  P1:=[op(P1),P[i]-1]:
fi:
od:

[S,[a+nops(S)-1,P1]]:
end:


#fPi(sPP): inputs a pair sPP=[S,[a1,P1]] where S is a set of integers, subset of [1,nops(P1)+nops(S)] and outputs
#a pair [a,P] where P is an a-parking function. It is the inverse of fP(P). Try:
#fPi(fP([1,[1,2,3,4]]));
fPi:=proc(SPP) local S, a1,P1,n,a,P,i:
S:=SPP[1]:
a1:=SPP[2][1]:
P1:=SPP[2][2]:
n:=nops(S)+nops(P1):
a:=a1-nops(S)+1:

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

if P1<>[] then
 RETURN(FAIL):
fi:
[a,P]:
end:


 


#gP(F): inputs a forest F (in new format [Emps,Bos,SetOfEdges]) and outputs the pair [S,F1] where
#S is the subset of Emps connected to the smallest member of Bos,let's call it b1, and F1 is the tree 
#[Emps minus S,Bos minus {b1} union S,SetOfEdges without b1]. Try:
#gP([{2,3,4},{1},{{1,2},{1,3},{1,4}}]);
gP:=proc(F) local S,b1,Emps,Bos,G,edg,sakh:
 Emps:=F[1]:
 Bos:=F[2]:
 G:=F[3]:
  b1:=min(op(Bos)):
 
  Bos:=Bos minus {b1}:

  S:={}:

 for edg in G do
  if member(b1,edg) then 
   G:=G minus {edg}:
   sakh:=edg minus {b1}:
   sakh:=sakh[1]:
    S:=S union {sakh}:
   Emps:=Emps minus  {sakh}:
   Bos:=Bos union {sakh}:
 fi:
od:

[S,[Emps,Bos,G]]:

end:

  
  
#gPi(Emps,Bos,P): inputs a set of Employees,Emps, a set of Bosses Bos, and a forest of rooted trees in new format, outputs the inverse of gP. Try
#gPi({2,3,4},{1}, gP([{2,3,4},{1},{{1,2},{1,3},{1,4}}])); 
gPi:=proc(Emps,Bos,F) local Emps1,Bos1,S,G1,b1,s:
 S:=F[1]:
 Emps1:=F[2][1]:
 Bos1:=F[2][2]:
 G1:=F[2][3]:

 if S minus Emps<>{} then
   RETURN(FAIL):
 fi:
 
  Emps1:=Emps1 union S:
  Bos1:=Bos1 minus S:
    if nops(Bos minus Bos1)<>1 then
      RETURN(FAIL):
    fi:
  b1:=(Bos minus Bos1)[1]:

  G1:=G1 union {seq({b1,s}, s in S)}:

[Emps,Bos,G1]:

end:


#CheckfP(n,a): checks that fPi is the inverse of fP (and vice versa), for all a-parking functions of length n. Try
#CheckfP(4,1);
CheckfP:=proc(n,a) local gu,mu,gu1,mu1:

gu:=PnaNew(n,a):

mu:={seq(fP(gu1),gu1 in gu)}:

evalb({seq(evalb(fPi(fP(gu1))=gu1),gu1 in gu)}={true}) and evalb({seq(evalb(fP(fPi(mu1))=mu1),mu1 in mu)}={true}):


end:





#CheckgP(Emps,Bos): checks that gPi is the inverse of gP (and vice versa), for all forests of rooted trees with set of roots Bos and
#set of NonRoots Emps. Try
#CheckgP({1,2},{3,4});
CheckgP:=proc(Emps,Bos) local gu,mu,gu1,mu1:

gu:=TnaEB(Emps,Bos):
mu:={seq(gP(gu1),gu1 in gu)}:

evalb({seq(evalb(gPi(Emps,Bos,gP(gu1))=gu1),gu1 in gu)}={true}) and evalb({seq(evalb(gP(gPi(Emps,Bos,mu1))=mu1),mu1 in mu)}={true}) :



end:



#SuS(S0,U): inputs a subset S0 of {1, ..., nops(U)} and outputs the corresponding subset of U. Try
#SuS({1,3},{4,5,7,8});
SuS:=proc(S0,U) local n,U1,i:
n:=nops(U):

if S0 minus {seq(i,i=1..n)}<>{} then
 RETURN(FAIL):
fi:

U1:=sort(convert(U,list)):
{seq(U1[S0[i]],i=1..nops(S0))}:
end:


#PtoTEB(P,Emps,Bos): inputs an a:=nops(Bos)-parking function and outputs the corresponding forest of rooted trees with a components.
#with a prescribed set of Employers and Bosses, Emps and Bos resepctively
#Try: 
#PtoTEB([2,[2,2,4]],{3,4,5},{1,2});
PtoTEB:=proc(P,Emps,Bos) local a,gu,S,lu,Bos1,Emps1,P1,n:
a:=P[1]:


if nops(Bos)<>a then
 RETURN(FAIL):
fi:

n:=nops(P[2]):


 if n<>nops(Emps) then
  RETURN(FAIL):
 fi:

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


gu:=fP(P):
S:=convert(gu[1],set):
S:=SuS(S,Emps):
P1:=gu[2]:

Emps1:=Emps minus S:
Bos1:=Bos minus {min(op(Bos))} union S:



lu:=PtoTEB(P1,Emps1,Bos1):


gPi(Emps,Bos,[S, lu]) :




end:


#CheckPtoTEB(n,a): Checks PtoT for a set of Employees with n elements and a set of bosses with a elements. Try:
#CheckPtoTEB(4,1);
CheckPtoTEB:=proc(n,a) local Emps, Bos,i,gu,P:
gu:=PnaNew(n,a):
Emps:={seq(i,i=1..n)}:
Bos:={seq(i,i=n+1..n+a)}:
evalb({seq(PtoTEB(P,Emps,Bos),P in gu)}=TnaEB(Emps,Bos)):
end:


#Milon(n): all the images of the 1-parking functions of length n as rooted trees with 0 the root. Try:
#Milon(4);
Milon:=proc(n) local gu,Emps,Bos,gu1,i,mu1:
gu:=PnaNew(n,1):
Bos:={0}:
Emps:={seq(i,i=1..n)}:

for gu1 in gu do
 mu1:=PtoTEB(gu1,Emps,Bos)[3]:
 print(gu1[2], mu1):
od:

end:

#MilonV(n): all the images of the 1-parking functions of length n as rooted trees with 0 the root, arranged by sorting, verbose version. Try:
#MilonV(3);
MilonV:=proc(n) local gu,mu,Bos,Emps,i,gu1,mu1:
Bos:={0}:
Emps:={seq(i,i=1..n)}:

gu:=Cna1(n,1):

for gu1 in gu do

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

print(`For the  permutations of `, gu1, ` we have `):

 print(gu1, PtoTEB([1,gu1],Emps,Bos)[3]):

 mu:=convert(permute(gu1),set) minus {gu1}:

 for mu1 in mu do

  print(mu1, PtoTEB([1,mu1],Emps,Bos)[3]):
 
 od:

od:

end:


#MilonVV(n): all the images of the 1-parking functions of length n as rooted trees with 0 the root, arranged by sorting, 
#verbose version, also draws them. Try:
#MilonVV(3);
MilonVV:=proc(n) local gu,mu,Bos,Emps,i,gu1,mu1:
Bos:={0}:
Emps:={seq(i,i=1..n)}:

gu:=Cna1(n,1):

for gu1 in gu do

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

print(`For the  permutations of `, gu1, ` we have `):

 print(gu1, PtoTEB([1,gu1],Emps,Bos)[3]):
 print(``):
print( DrawRT( PtoTEB([1,gu1],Emps,Bos)[3],0));

 mu:=convert(permute(gu1),set) minus {gu1}:

 for mu1 in mu do

  print(mu1, PtoTEB([1,mu1],Emps,Bos)[3]):
 print(``):
  print(DrawRT( PtoTEB([1,mu1],Emps,Bos)[3],0));
 print(``):
 
 od:

od:

end:



#Roul(L): inputs a list L and outputs i from 1 to nops(L) with prob. L[i]/convert(L,`+`); Try
#Roul([1,2,1]);
Roul:=proc(L) local i,su,r,i1:
su:=convert(L,`+`):

r:=rand(1..su)():

for i from 1 while add(L[i1],i1=1..i-1)<r do od:
i-1:
end:


#RW(a,b): a random walk from [0,0] to [a,b] that ends at [a,b]. Try:
#RW(2,4);
RW:=proc(a,b) local gu,ro:
if a=0 and b=0 then
 RETURN([]):
fi:

ro:=Roul([a,b]):
if ro=1 then
 gu:=RW(a-1,b):
 RETURN([op(gu),0]):
else
  gu:=RW(a,b-1):
 RETURN([op(gu),1]):
fi:
end:




#RSU(n,k): a random k-subset  of {1, ..., n}. Try: RSU(5,3);
RSU:=proc(n,k) local gu,mu,i:
gu:=RW(n-k,k):

mu:={}:

for i from 1 to n do
 if gu[i]=1 then
   mu:=mu union {i}:
 fi:
od:

mu:

end:


#RPF(n,a): a random a-parking function of length n. Try
#RPF(5,1);
RPF:=proc(n,a) local k,P1,P,S,i:

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

if n=1 and a=1 then
 RETURN([1]):
fi:

k:=Roul([seq(binomial(n,k)*tn(n-k,a+k-1),k=0..n)])-1:


S:=RSU(n,k):


P1:=RPF(n-k,a+k-1):


P:=[]:
for i from 1 to n do

 if member(i,S) then
  P:=[op(P),1]:
 else
  P:=[op(P),P1[1]+1]:
  P1:=[op(2..nops(P1),P1)]:
 fi:

od:

P:

end:

  
#PtoT0(P): inputs a 1-parking function and outputs the rooted tree with root 0 and non-roots labelled {1, ..., nops(P)}. 
PtoT0:=proc(P) local i:
PtoTEB([1,P],{seq(i,i=1..nops(P))},{0})[3] :
end:

#PtoT1(P): inputs a 1-parking function and outputs the rooted tree with root 1 and non-roots labelled {2, ..., nops(P)+1}. 
PtoT1:=proc(P) local i:
PtoTEB([1,P],{seq(i,i=2..nops(P)+1)},{1})[3] :
end:

#RRT(n): a random labelled rooted tree with root 0 and non-roots labelled {1, ..., n}. Try:
#RRT(10);
RRT:=proc(n) :
PtoT0(RPF(n,1)):
end:


#Neighs(G,v): given a graph G as a set of edges, finds all the neighbors of v, followed by the set of edges. Try:
#Neighs({{1,2},{1,3},{3,4},{1,5}},1);
Neighs:=proc(G,v) local edg,gu,mu:

gu:={}:
mu:={}:
for edg in G do
 if member(v,edg) then
  mu:=mu union {edg}:
   gu:=gu union {(edg minus {v})[1]}:
 fi:
od:

gu,mu:

end:

#CC(G,v): the connected component of the vertex v in the graph G. Try:
#CC({{1,2},{1,3},{3,4},{1,5}},1);
CC:=proc(G,v) local gu,gu1,gu2,v1,mu:
gu:={v}:
gu1:=gu union Neighs(G,v)[1]:
mu:={}:
while gu<>gu1 do
  mu:=mu union {seq(op(Neighs(G,v1)[2]),v1 in  gu1)}:
 gu2:=gu1 union {seq(op(Neighs(G,v1)[1]),v1 in  gu1)}:
 gu:=gu1:
 gu1:=gu2:
od:
gu1,mu:

end:


#CanF(T,r): the canonical form of a rooted tree T with root r. Try
#CanF({{0,1},{0,2},{0,3},{1,4}},0);
CanF:=proc(T,r) local gu,yel,T1,i,lu:
gu:=CC(T,r);

if gu[1]={r} then
 RETURN([r,[]]):
fi:

yel:=Neighs(T,r):

T1:=T minus yel[2]:

yel:=yel[1]:
yel:=sort(convert(yel,list)):

gu:=[]:

 for i from 1 to nops(yel) do
  lu:=CC(T1,yel[i])[2]:
  lu:=CanF(lu,yel[i]):
  gu:=[op(gu),lu]:
 od:

[r,gu]:

end:

#Dorot(T): inputs a rooted tree, and outputs, in order, the successive generations. Try:
#Dorot(CanF({{0,1},{0,2},{0,3},{1,4}},0));
Dorot:=proc(T) local gu,mu,lu,i,omek,L,j,ku:


if T[2]=[] then
 RETURN([[T[1]]]):
fi:
gu:=[[T[1]]]:

lu:=[seq(T[2][i][1],i=1..nops(T[2]))]:
gu:=[op(gu),lu]:

for i from 1 to nops(T[2]) do
 ku:=Dorot(T[2][i]):
 ku:=[op(2..nops(ku),ku)]:
 L[T[2][i]]:=ku:
od:

omek:=max(seq(nops(L[T[2][i]]),i=1..nops(T[2]))):

for j from 1 to omek do
 mu:=[]:
  for i from 1 to nops(T[2]) do
    if  nops(L[T[2][i]])>=j then
       mu:=[op(mu),op(L[T[2][i]][j])]:
   fi:
  od:
 gu:=[op(gu),mu]:
od:
 
gu: 
end:


#DrawRT(T,r): draws the rooted tree T with root r. Try
#DrawRT({{0,1},{0,2},{0,3}},0);
DrawRT:=proc(T,r) local S,lu,gu,d,j,i,k,eps:
eps:=0.3:
S:=CanF(T,r):
lu:=Dorot(S):

if nops(lu[1])<>1 then
  RETURN(FAIL):
fi:

d[1]:=-1/2:

    gu:=plot([[d[1],0]],style=point,axes=none,scaling=constrained,view=[-3..3,-5..1]):
 
    gu:=gu,textplot([d[1]+eps,0,lu[1][1]]):

for i from 2 to nops(lu) do
d[i]:=-nops(lu[i])/2:
 for j from 1 to nops(lu[i]) do

    gu:=gu,plot([[d[i]+j-1,-i+1]],style=point):

   gu:=gu,textplot([d[i]+j-1+eps,-i+1, lu[i][j]]):

   for k from 1 to nops(lu[i-1]) do
     if member({lu[i][j],lu[i-1][k]},T) then
       gu:=gu,plot([[d[i]+j-1,-i+1], [d[i-1]+k-1,-i+2]]):
     fi:
   od:
 od:
od:

display(gu);

end: