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


read `YBMdata.txt`:

print(`First Written: June 2022: tested for Maple 2020 `):
print(`Version  : June 2022  `):
print():
print(`This is RedDec.txt, A Maple package`):
print(`accompanying Shalsoh B. Ekhad and Doron Zeilberger's  article: `):
print(` Experimenting with reduced decompositions in the Symmetric Group`):

print():
print(`The most current version is available on WWW at:`):
print(` http://sites.math.rutgers.edu/~zeilberg/tokhniot/RedDec.txt .`):
print(`Please report all bugs to: DoronZeil at gmail dot com .`):

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(`For a list of the Edeleman-Greene functions type: ezraEG();`):
print():
 




ezraEG:=proc()
if args=NULL then

print(`The Edeleman-Greene procedures are`):
print(`  EG, EGr, jeuT, MilimS, RandW, SchP, TtoW, WtoT,  `):

else
ezra(args):

fi:


end:

ezra1:=proc()
if args=NULL then

print(`The SUPPORTING procedures are`):
print(` AllYF, Alpha, Conj, CycS1, CycS, DecXnij, Des, EvalMila, EXnD,`):
print(`GNW, inv, InvTab, jeuT, KickOut, LocatePlaces,MilaToCh, Mul, Pars, Pars1, ParsL, PermuteList, PrT, PrXnij,  PrXnrj1j2, RandCh,Shn, Shnj, Sh1, Sh2, tin, w0n, w0nj, YB, YBp, YF `):

else
ezra(args):

fi:


end:

ezra:=proc()
if args=NULL then

print(` RedDec.txt: A Maple package for  `):

print(`The MAIN procedures are:  EXn, GFm, Gp, GpA, GNWs, Milim, NuMilim, Nus, NuXnij, RandCh1, ReinerTable, Simu, SYT, Xnij, Xni1j1i2j2, Xnrj1j2, YBP `):
print(``):



elif nargs=1 and args[1]=AllYF then
print(`AllYF(n): All the YF(L) for L in the stair-case shape [n-1,n-2,...,1]. Try:`):
print(`AllYF(4);`):

elif nargs=1 and args[1]=Alpha then
print(`Alpha(f,x,N): Given a probability generating function`):
print(`outputs statistical information. Try:`):
print(`Alpha(((1+x)/2)^100,x,4);`):

elif nargs=1 and args[1]=Conj then
print(`Conj(L): The conjugate of the partition L. Try:`):
print(`Conj([5]);`):


elif nargs=1 and args[1]=CycS then
print(`CycS(w,n,r): The cyclic shift to the left of the word w, by r paces, in the alphabet {1,...,n-1}. Try: `):
print(`CycS([1,2,1],3,2);`):

elif nargs=1 and args[1]=CycS1 then
print(`CycS1(w,n): The cyclic shift of the word w in the alphabet {1,...,n-1}. Try: `):
print(`CycS1([1,2,1],3);`):

elif nargs=1 and args[1]=DecXnij then
print(`DecXnij(n,i,j): The set of pairs [pi1,pi2] such that inv(pi1)=i-1, inv(pi2)=n*(n-1)/2-i-2 and [n,...,1]=pi1 [j+2,j+1,j] pi2, via Xnij(n,i,j). Try:`):
print(`DecXnij(5,1,1);`):

elif nargs=1 and args[1]=Des then
print(`Des(pi): The set of descentes of the permutation pi. Try:`):
print(`Des([5,3,1,4,2]);`):

elif nargs=1 and args[1]=EG then
print(`EG(w,n):  Fiven a word w  in {1,...,n-1} of length n*(n-1)/2 that represents the longest permutation [n,...,1] applies the Edelman-Greene mapping and outputs a balanced`):
print(`tableau of shaple [n-1,n-2,...2,1]. Try:`):
print(`w:=RandCh1(10)[2]: EG(w,10);`):


elif nargs=1 and args[1]=EGr then
print(`EGr(T,n):  Fiven a balanced tableau of stairase shape [n-1,...,1], applies the reverse Edelman-Greene mapping and outputs a word in the alphabet {1, ..., n-1}`):
print(`of length n(n-1)/2 describing a chain in the weak Bruhat order. try:`):
print(`w:=RandCh1(10)[2]: EGr(EG(w,10),10);`):


elif nargs=1 and args[1]=EvalMila then
print(`EvalMila(w,n): Given a word w=[w1, ..., wk], in {1,...,n-1} finds the product s_w[1].... s_w[k] Try: `):
print(`EvalMila([1,3,2,3],4);`):

elif nargs=1 and args[1]=EXn then
print(`EXn(n): The expected number of Yang-Baxter moves on the symmetric group on n elements. Try:`):
print(`EXn(8);`):


elif nargs=1 and args[1]=EXnD then
print(`EXnD(n): Like EXn(n) but done directly. Only for checking. Try:`):
print(`EXnD(6);`):

elif nargs=1 and args[1]=GFm then
print(`GFm(pi,z1,z2): The generating function of the reduced words that evaluate to pi according to the weight z1^YB(w)[1]*z2^YB(w)[2]`):

elif nargs=1 and args[1]=GNW then
print(`GNW(L): A random Standard Young tableau of shape L`):
print(`using the Greene-Nijenhus-Wilf algorithm. For example do:`):
print(`GNW([3,3,3]);`):


elif nargs=1 and args[1]=GNWs then
print(`GNWs(n): A random Standard Young tableau of shape [n-1,...,1]`):
print(`using the Greene-Nijenhus-Wilf algorithm. For example do:`):
print(` GNWs(10); `):

elif nargs=1 and args[1]=Gp then
print(`Gp(S): The group generated by the list of permutations S. Try:`):
print(` Gp([[2,1,3],[1,3,2]]);`):

elif nargs=1 and args[1]=GpA then
print(`GpA(S,a): All the members  of Gp(S) whose minimal word has length A. Try:`):
print(`GpA([[2,1,3],[1,3,2]],2);`):

elif nargs=1 and args[1]=inv then
print(`inv(pi): The number of inversions of the permutation pi. Try:`):
print(`inv([5,1,3,2,4]);`):

elif nargs=1 and args[1]=InvTab then
print(`InvTab(pi): The inversion Table  of the permutation pi. Try:`):
print(`InvTab([5,1,2,3,4]);`):

elif nargs=1 and args[1]=jeuT then
print(`jeuT(T,n,k): Applying jeu-de-taquin to a standard tableu T of k+1,...,k+n . Try:`):
print(`jeuT([[1,3,4],[2,6],[5]],6,0);`):

elif nargs=1 and args[1]=KickOut then
print(`KickOut(w,i,j,n): Inputs a word w in {1,...,n-1}, kicks out w[i]...w[j] and outputs the resulting permutation`):
print(`and its Sh1, and Sh2. try:`):
print(`KickOut([1, 2, 1, 4, 3, 2, 4, 1, 3, 2],1,3,5);`):

elif nargs=1 and args[1]=LocatePlaces then
print(`LocatPlaces(T,S): Given a tableau T and a set of entries S, outputs the set of places where they reside. Try:`):
print(`LocatePlaces([[1,2,3,4],[5,6,7],[8,9]],{1,6});`):

elif nargs=1 and args[1]=MilaToCh then
print(`MilaToCh(W,n): converts the word W in S_n to a chain. Try:`):
print(`MilaToCh([1,2,1],3);`):

elif nargs=1 and args[1]=Milim then
print(`Milim(pi): All the minimal ways of expressing pi as a prodcut of trasnpositions of the type [i+1,i]`):

elif nargs=1 and args[1]=MilimS then
print(`MilimS(n): all the words representing minimal factorizations of [n,n-1,..,2,1], using TtoW(T,n). Try:`):
print(`MilimS(5);`):

elif nargs=1 and args[1]=Mul then
print(`Mul(a,b): The product of permutation a and b. Try:`):
print(`Mul([3,1,2,4],[4,2,1,3]);`):


elif nargs=1 and args[1]=NuMilim then
print(`NuMilim(pi): the cardinality of Milim(pi) (q.v.). Try:`):
print(`NuMilim([4,3,2,1]);`):

elif nargs=1 and args[1]=Nus then
print(`Nus(n): The cardinalities that show up in the cardinalities of the intersections of Xnij(n1,i1,j1) and Xnij(n2,i2,j2). Try:`):
print(`Nus(4);`):

elif nargs=1 and args[1]=NuXnij then
print(`NuXnij(n,i,j): Same as nops(Xnij(n,i,j)) but done cleverly. Try:`):
print(`NuXnij(5,1,1);`):

elif nargs=1 and args[1]=Pars then
print(`Pars(n,k): the set of partitions of . Try:`):
print(`Pars(5);`):


elif nargs=1 and args[1]=ParsL then
print(`ParsL(L): the set of sub-shapes of L . Try:`):
print(`ParsL([3,2,1]);`):

elif nargs=1 and args[1]=Pars1 then
print(`Pars1(n,k): the set of partitions of n with largest part k. Try:`):
print(`Pars1(5,2):`):

elif nargs=1 and args[1]=PermuteList then
print(`PermuteList(n): The list of length binomial(n,2)+1 such that for i=0..binomial(n,2), L[i+1] is the set of permutations of length n with i inversions. Try:`):
print(`PermuteList(5);`):

elif nargs=1 and args[1]=PrT then
print(`PrT(T): prints the tableau T. Try:`):
print(`PrT([[3,2,1],[4,5],[6]]);`):

elif nargs=1 and args[1]=PrXnij then
print(`PrXnij(n,i,j): Prints out the Tableaux images under the Edelman-Greene mapping of the words given by Xnij(n,i,j) (q.v.). Try:`):
print(`PrXnij(4,1,1);`):


elif nargs=1 and args[1]=PrXnrj1j2 then
print(`PrXnrj1j2(n,r,j1,j2): Prints out the Tableaux images under the Edelman-Greene mapping of the words given by Xnrj1j2(n,r,j1,j2) (q.v.) Try:`):
print(`PrXnrj1j2(5,4,1,2);`):

elif nargs=1 and args[1]=RandCh then
print(`RandCh(pi): A random maximal chain in the weak Bruhat order from the permutation pi to the identity. It is not "uniformaly at random". It is just to generate`):
print(`many examples. It also outputs the correponding reduced decomposition of pi. Try:`):
print(`RandCh([5,4,3,2,1]);`):


elif nargs=1 and args[1]=RandCh1 then
print(`RandCh1(n): A random maximal chain in the weak Bruhat order from the permutation [n,...,1] to the identity. It is not "uniformaly at random". It is just to generate examples.`):
print(`It returns the actual chain with binomial(n,2)+1 links, and the implied word of length binomial(n,2).`):
print(`Try: `):
print(`RandCh1(20);`):

elif nargs=1 and args[1]=RandW then
print(`RandW(n): a uniformly at random word representing [n,n-1,...,1]. Try:`):
print(`RandW(10);`):

elif nargs=1 and args[1]=ReinerTable then
print(`ReinerTable(w): Inputs a permutation w (of length n=nops(w)) and outputs the list 0f lists of length n-2 where each list has length inv(w)-2 , let's call it L such that`):
print(`L[j][k] is the set of minimal words that evaluate to w that have either [j,j+1,j] or [j+1,j,j+1] starting at the k-th place. Try:`):
print(`ReinerTable([5,4,3,2,1]);`):


elif nargs=1 and args[1]=SchP then
print(`SchP(T,n)) applies the Schutzenberger transform to tableu T of {1,...,n}. Try:`):
print(`SchP([[1,3,4],[2,6],[5]]);`):

elif nargs=1 and args[1]=Sh1 then
print(`Sh1(w): What Richard Stanley calls Lambda(w) in Eur. J. Combinatorics 5 (1984), 359-372, p. 367. Try:`):
print(`Sh1([4,1,3,2]);`):

elif nargs=1 and args[1]=Sh2 then
print(`Sh2(w): What Richard Stanley calls mu(w) in Eur. J. Combinatorics 5 (1984), 359-372, p. 367. Try:`):
print(`Sh2([4,1,3,2]);`):

elif nargs=1 and args[1]=Shn then
print(`Shn(n): the partition [n-1,...,1]. Try:`):
print(`Shn(8);`):

elif nargs=1 and args[1]=Shnj then
print(`Shnj(n,j): the partition [n-1,...,1]- [0^(j-1),2,1,0^(n-2-j)]. It is the shape of w0nj(n,j).  Try:`):
print(`Shnj(10,1);`):

elif nargs=1 and args[1]=Simu then
print(`Simu(n,K,L): generates K random maximal words  representing [n,n-1,...,1] and retruns the list whose i-th item is the ratio of those that have i-1 Yang-Baxter moves (up to the max), followed`):
print(`by a list of length L whose first item is the average, second the variance, followed by the moments about the mean. Try:`):
print(`Simu(10,100,4);`):

elif nargs=1 and args[1]=SYT then
print(`SYT(L): all the standard Young Tableaux of shape L`):
print(`For example, try: SYT([2,2]);`):

elif nargs=1 and args[1]=tin then
print(`tin(i,n): the transposition (i,i+1) of {1,...,n}. Try: tin(6,2);`):

elif nargs=1 and args[1]=TtoW then
print(`TtoW(T,n): inputs a Young tableau of staircase shape [n-1,...,1] outputs a word in {1,..,n-1} that is a minimal decomposition of the permutation [n,n-1,..,1]. Try:`):
print(`TtoW([[1,2,3,4],[5,6,7],[8,9],[10]],5);`):

elif nargs=1 and args[1]=w0n then
print(`w0n(n,j): The permutation [n,n-1,n-2,..,1]. Try:`):
print(`w0n(10);`):

elif nargs=1 and args[1]=w0nj then
print(`w0nj(n,j): The permutation (j+2,j+1,j)*[n,n-1,...,1] in other words [n,n-1,..., j+3,j,j+1,j+2,j-1,...,1]. Try:`):
print(`w0nj(10,2);`):

elif nargs=1 and args[1]=WtoT then
print(`WtoT(w,n): inputs a tableau T of staircase  shape [n-1,...,1] and outputs the corresponding word representing a minimal factorization of the permutation [n,n-1,...,1]. According to`):
print(`Edelman-Greene. Try:`):
print(`WtoT([1,2,1],3);`):

elif nargs=1 and args[1]=Xnij then
print(`Xnij(n,i,j): Inputs a positive integer n, and positive integer i between 1 and binomial(n,2)-2 and j such that 1<=j<=n-2, outputs The set of reduced words w of [n,...,1] that have the property `):
print(`w[i]=j,w[i+1]=j+1, w[i+2]=j or w[i]=j+1,w[i+1]=j, w[i+2]=j+1 . Try:`):
print(`Xnij(5,1,1);`):


elif nargs=1 and args[1]=Xnrj1j2 then
print(`Xnrj1j2(n,r,j1,j2): The set of minimal words of w0n(n) with Yang-Baxter moves at location 1 of type j1, and at location r of type j2. Try:`):
print(`Xnrj1j2(6,3,1,2);`):

elif nargs=1 and args[1]=Xni1j1i2j2 then
print(`Xni1j1i2j2(n,i1,j1,i2,j2): The set of minimal words of w0n(n) with Yang-Baxter moves at location i1 of type j1, and at location i2 of type j2. Try:`):
print(`Xni1j1i2j2(6,1,1,9,1);`):

elif nargs=1 and args[1]=YB then
print(`YB(w): The number of Yang-Baxter moves in a word w of {1,...,n-1}. It returns the pair [co1,co2] where co1 is the type [j,j+1,j] and co2 the type [j+1,j,j+1].`):


elif nargs=1 and args[1]=YBp then
print(`YBp(w): The sets where Yang-Baxter moves in a word w of {1,...,n-1} occur. It returns the pair [co1,co2] where co1 is the set of places that start [j,j+1,j] and co2 for [j+1,j,j+1].`):

elif nargs=1 and args[1]=YBP then
print(`YBP(n,z1,z2): The generating function for the longest word. In other words:`):
print(`GFm([seq(n+1-i,i=1..n)],z1,z2); Try:`):
print(`YBP(5,z,z);`):

elif nargs=1 and args[1]=YF then
print(`YF(L): The Young-Frobenius formula. Try:`):
print(`YF([5,3,2]);`):

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

fi:


end:


with(combinat):


###START FROM GreeneNiejenhuisWilf

#OneStepGNW1(L,cell): Given a shape L and a cell=[i,j]
#decides where in the hook to go
OneStepGNW1:=proc(L,cell) local H,i,j,i1,j1:

i:=cell[1]:
j:=cell[2]:

if not (1<=i and i<=nops(L) and 1<=j and j<=L[i]) then
 ERROR(`Bad input`):
fi:

H:=[seq([i,j1],j1=j+1..L[i])]:

for i1 from i+1 to nops(L) while j<=L[i1] do
 H:=[op(H),[i1,j]]:
od:

if H=[] then
 RETURN(FAIL):
fi:

H[rand(1..nops(H))()]:
end:




#OneStepGNW(L): Given a shape L applies one step
#of the Greene-Nijenhius-Wilf algorithm to decide
#where to put n. Try:
#OneStepGNW([2,2]);
OneStepGNW:=proc(L) local cell,n,cell1,cell2,T1,i1,j1:
n:=convert(L,`+`):
T1:=[seq(seq([i1,j1],j1=1..L[i1]),i1=1..nops(L))]:
if n=1 then
 RETURN([1,1]):
fi:


cell:=T1[rand(1..n)()]:

cell1:=OneStepGNW1(L,cell):

if cell1=FAIL then
 RETURN(cell):
fi:

while cell1<>FAIL do
cell2:=OneStepGNW1(L,cell1):
cell:=cell1:
cell1:=cell2:
od:
cell:

end:


#GNW(L): A random Standard Young tableau of shape L
#using the Greene-Nijenhus-Wilf algorithm. For example do:
#GNW([3,3,3]);
GNW:=proc(L) local L1,cell,T1,n,i:

n:=convert(L,`+`):
if n=1 then
RETURN([[1]]):
fi:

cell:=OneStepGNW(L):
i:=cell[1]:

if L[i]>1 then
L1:=[op(1..i-1,L),L[i]-1,op(i+1..nops(L),L)]:
else
L1:=[op(1..i-1,L),op(i+1..nops(L),L)]:
fi:


T1:=GNW(L1):

if L[i]>1 then
[op(1..i-1,T1),[op(T1[i]),n],op(i+1..nops(T1),T1)]:
else
[op(1..i-1,T1),[n]]:
fi:

end:

#SYT(L): all the standard Young Tableaux of shape L
#For example, try: SYT([2,2]);
SYT:=proc(L) local gu,n,L1,mu,mu1,i:
option remember:

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

n:=convert(L,`+`):

gu:={}:

for i from 1 to nops(L) do

if i=nops(L) or L[i]>L[i+1] then
 L1:=[op(1..i-1,L),L[i]-1,op(i+1..nops(L),L)]:
 if L1[nops(L1)]=0 then
  L1:=[op(1..nops(L1)-1,L1)]:
 fi:
mu:=SYT(L1):

if nops(L1)<nops(L) then
 gu:=gu union {seq([op(mu1),[n]],mu1 in mu)}:
else
 gu:=gu union 
{
seq(
[op(1..i-1,mu1),
[op(mu1[i]),n],
op(i+1..nops(mu1),mu1)],
mu1 in mu)
}:
fi:

fi:
od:

gu:

end:


###END FROM GreeneNiejenhuisWilf

#Start from Histabrut
#AveAndMoms(f,x,N): Given a probability generating function
#f(x) (a polynomial, rational function and possibly beyond)
#returns a list whose first entry is the average 
#(alias expectation, alias mean)
#followed by the variance, the third moment (about the mean) and
#so on, until the N-th moment (about the mean).
#If f(1) is not 1, than it first normalizes it by dividing
#by f(1) (if it is not 0) .
#For example, try:
#AveAndMoms(((1+x)/2)^100,x,4);

AveAndMoms:=proc(f,x,N) local mu,F,memu1,gu,i:
mu:=simplify(subs(x=1,f)):

if mu=0 then
print(f, `is neither a prob. generating function nor can it be made so`):
RETURN(FAIL):
fi:

F:=f/mu:


memu1:=simplify(subs(x=1,diff(F,x))):

gu:=[memu1]:

F:=F/x^memu1:

F:=x*diff(F,x):
for i from 2 to N do
 F:=x*diff(F,x):
 gu:=[op(gu),simplify(subs(x=1,F))]:
od:

gu:

end:


#Alpha(f,x,N): Given a probability generating function
#f(x) (a polynomial, rational function and possibly beyond)
#returns a list, of length N, whose 
#(i) First entry is the average
#(ii): Second entry is the variance
#for i=3...N, the i-th entry is the so-called alpha-coefficients
#that is the i-th moment about the mean divided by the
#variance to the power i/2 (For example, i=3 is the skewness
#and i=4 is the Kurtosis)
#If f(1) is not 1, than it first normalizes it by dividing
#by f(1) (if it is not 0) .
#For example, try:
#Alpha(((1+x)/2)^100,x,4);
Alpha:=proc(f,x,N) local gu,i:
 gu:=AveAndMoms(f,x,N):
 if gu=FAIL then
  RETURN(gu):
 fi:

if gu[2]=0 then
print(`The variance is 0`):
  RETURN(FAIL):
fi:

[gu[1],gu[2],seq(gu[i]/gu[2]^(i/2),i=3..N)]:

end:

#end from Histabrut

#Mul(a,b): The product of permutation a and b
Mul:=proc(a,b) local i:
[seq(b[a[i]],i=1..nops(a))]:
end:

#Gp(S): The group generated by the list of permutations S. Try: Gp([[2,1,3],[1,3,2]]):
Gp:=proc(S) local G,n,i,N,N1,s,g:
option remember:
n:=nops(S[1]):

G:={[seq(i,i=1..n)]}:

N:=G union {seq(seq(Mul(s,g),g in G), s in S)}:

while G<>N do

N1:=N union {seq(seq(Mul(s,g),g in N), s in S)}:

G:=N:
N:=N1:
od:
G:
end:


#GpA(S,a): All the members  of Gp(S) whose minimal word has length A. Try:
#GpA([[2,1,3],[1,3,2]],2);
GpA:=proc(S,A) local n,i,katan,g,muam,mu,gu:
option remember:

n:=nops(S[1]):

if A=0 then
 RETURN({[seq(i,i=1..n)]}):
fi:

gu:=GpA(S,A-1):
katan:={seq(op(GpA(S,i)),i=0..A-1)}:

mu:={}:
for g in gu do
 for i from 1 to nops(S) do
  muam:=Mul(g,S[i]):
  if not member(muam,katan) then
    mu:=mu union {muam}:
  fi:
 od:
od:

mu:

end:


#inv(pi): The number of inversions of the permutation pi
inv:=proc(pi) local n,co,i,j:
n:=nops(pi):
co:=0:


for i from 1 to n do
 for j from i+1 to n do
 if pi[i]>pi[j] then
   co:=co+1:
 fi:
 od:
od:
co:
end:




#Sh1(w): What Stanley calls Lambda(w) in Europ. J. Combinatorics 5 (1984), 359-372 The conjugate of the sorted inversion table 
Sh1:=proc(w) local n,i,j,gu,gu1:
n:=nops(w):


gu:=[]:
for i from 2 to n do
gu1:=0:
 for j from 1 to i-1 do
 if w[j]>w[i] then
   gu1:=gu1+1:
 fi:
 od:
gu:=[op(gu),gu1]:
od:

gu:=sort(gu,`>`):

for i from 1 to nops(gu) while gu[i]<>0 do od:
[op(1..i-1,gu)]:
end:

#Sh2(w): What Stanley calls Mu(w) in Europ. J. Combinatorics 5 (1984), 359-372 The conjugate of the sorted inversion table 
Sh2:=proc(w) local n,i,j,gu,gu1:
n:=nops(w):


gu:=[]:
for i from 1 to n-1 do
gu1:=0:
 for j from i+1 to n do
 if w[j]<w[i] then
   gu1:=gu1+1:
 fi:
 od:
gu:=[op(gu),gu1]:
od:
gu:=sort(gu,`>`):

for i from 1 to nops(gu) while gu[i]<>0 do od:
Conj([op(1..i-1,gu)]):
end:



#MilimOld(pi): All the minimal ways of expressing pi as a prodcut of trasnpositions of the type [i+1,i]
MilimOld:=proc(pi) local n,i,pi1,mu,lu,w:
option remember:

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

n:=nops(pi):

mu:={}:

for  i from 1 to n-1 do
 pi1:=[op(1..i-1,pi),pi[i+1],pi[i],op(i+2..n,pi)]:
  if inv(pi1)<inv(pi) then
   lu:=MilimOld(pi1):
    mu:=mu union {seq([i,op(w)],w in lu)}:
  fi:
od:
mu:


end:



#Milim(pi): All the minimal ways of expressing pi as a prodcut of trasnpositions of the type [i+1,i]
Milim:=proc(pi) local n,i,pi1,mu,lu,w:
option remember:

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

n:=nops(pi):

mu:={}:

for  i from 1 to n-1 do
if pi[i]>pi[i+1] then
 pi1:=[op(1..i-1,pi),pi[i+1],pi[i],op(i+2..n,pi)]:
   lu:=Milim(pi1):
    mu:=mu union {seq([i,op(w)],w in lu)}:
  fi:
od:
mu:


end:


#YB(w): The number of Yang-Baxter moves in a word w of {1,...,n-1}. It returns the pair [co1,co2] where co1 is the type [j,j+1,j] and co2 the type [j+1,j,j+1]
YB:=proc(w) local co1,co2,i:

co1:=0:
co2:=0: 
for i from 1 to nops(w)-2 do

  if w[i]=w[i+2] and w[i+1]-w[i]=1 then
   co1:=co1+1:
  elif w[i]=w[i+2] and w[i]-w[i+1]=1 then
   co2:=co2+1:
 fi:
od:


[co1,co2]:

end:

YBtot:=proc(w) local lu: lu:=YB(w): lu[1]+lu[2]:end:


#YBp(w): The sets of Yang-Baxter moves in a word w of {1,...,n-1}. It returns the pair [co1,co2] where co1 is the type [j,j+1,j] and co2 the type [j+1,j,j+1]
YBp:=proc(w) local co1,co2,i:

co1:={}:
co2:={}: 
for i from 1 to nops(w)-2 do

  if w[i]=w[i+2] and w[i+1]-w[i]=1 then
   co1:=co1 union {i}:
  elif w[i]=w[i+2] and w[i]-w[i+1]=1 then
   co2:=co2 union {i}:
 fi:
od:


[co1,co2]:

end:


#YBpp(w): The sets of Yang-Baxter moves in a word w of {1,...,n-1}. It returns the pair [co1,co2] where co1 is the type [j,j+1,j] and co2 the type [j+1,j,j+1]
YBpp:=proc(w) local co1,co2,i:

co1:={}:
co2:={}: 
for i from 1 to nops(w)-2 do

  if w[i]=w[i+2] and w[i+1]-w[i]=1 then
   co1:=co1 union {[i,w[i]]}:
  elif w[i]=w[i+2] and w[i]-w[i+1]=1 then
   co2:=co2 union {[i,w[i+1]]}:
 fi:
od:


[co1,co2]:

end:

#GFm(pi,z1,z2): The generating function of the reduced words that evaluate to pi according to the weight z1^YB(w)[1]*z2^YB(w)[2]
GFm:=proc(pi,z1,z2) local gu,w,co,lu:

gu:=Milim(pi):

co:=0:

for w in gu do
 lu:=YB(w):
 co:=co+z1^lu[1]*z2^lu[2]:
od:
co:
end:

#YBP(n,z1,z2): The generating function for the longest word. In other words:
#GFm([seq(n+1-i,i=1..n)],z1,z2); Try:
#YBP(5,z,z);
YBP:=proc(n,z1,z2) local i:
GFm([seq(n+1-i,i=1..n)],z1,z2):
end:



#YF(L): The Young-Frobenius formula
YF:=proc(L) local i,j,n,k:
n:=convert(L,`+`):
k:=nops(L):
n!/mul((L[i]+k-i)!,i=1..k)*mul(mul((L[j]+k-j)-(L[i]+k-i),j=1..i-1),i=1..k):
end:



#Des(pi): The set of descentes of the permutation pi
Des:=proc(pi) local S,i:
S:={}:

for i from 1 to nops(pi)-1 do
 if pi[i]>pi[i+1] then
  S:=S union {i}:
 fi:
od:

S:


end:

#RandCh(pi): A random maximal chain in the weak Bruhat order from the permutation pi to the identity. It is not "uniformaly at random". It is just to generate
#many examples. It also outputs the correponding reduced decomposition of pi. Try:
#RandCh([5,4,3,2,1]);
RandCh:=proc(pi) local S,i,C,n,i1,pi1,C1:
n:=nops(pi):


C:=[pi]:
C1:=[]:
pi1:=pi:

while pi1<>[seq(i1,i1=1..n)] do

 S:=Des(pi1):

i:=S[rand(1..nops(S))()]:

C1:=[op(C1),i]:
pi1:=[op(1..i-1,pi1),pi1[i+1],pi1[i],op(i+2..n,pi1)]:

C:=[op(C),pi1]:

od:

C,C1:
#[seq(C[nops(C)+1-i],i=1..nops(C))],[seq(C1[nops(C1)+1-i],i=1..nops(C1))]:

end:


#EvalMila(w,n): Given a word w=[w1, ..., wk], in {1,...,n-1} finds the product s_w[1].... s_w[k] Try
#EvalMila([1,3,2,3],4);
EvalMila:=proc(w,n) local gu,i,j,i1:
gu:=[seq(i,i=1..n)]:

for i from 1 to nops(w) do
 j:=w[i]:
 gu:=Mul(gu,[seq(i1,i1=1..j-1),j+1,j,seq(i1,i1=j+2..n)]):
od:
gu:
end:







RandCh1:=proc(n) local i:
RandCh([seq(n+1-i,i=1..n)]):
end:


#EG(w,n):  Fiven a word w  in {1,...,n-1} of length n*(n-1)/2 that represents the longest permutation [n,...,1] applies the Edelman-Greene mapping and outputs a balanced
#tableau of shaple [n-1,n-2,...2,1]. Try:
#w:=RandCh1(10)[2]: EG(w,n)
EG:=proc(w,n) local i,pi,k,b,r,j:

if nops(w)<>n*(n-1)/2 and EvalMila(w,n)<>[seq(n+1-i,i=1..n)] then
 RETURN(FAIL):
fi:

pi:=[seq(i,i=1..n)]:

for k from 1 to nops(w) do

 r:=w[k]:

 j:=pi[r]:
 i:=pi[r+1]:

 b[n+1-i,j]:=k:
 pi:=[op(1..r-1,pi),i,j,op(r+2..n,pi)]:
od:
[seq([seq(b[i,j],j=1..n-i)],i=1..n-1)]:
end:


#PrT(T): prints the tableau T
PrT:=proc(T) local i:
for i from 1 to nops(T) do
lprint(op(T[i])):
od:
end:


#Conj(L): The conjugate of the partition L. Try
#Conj([5]);

Conj:=proc(L) local L1,k,i1,i:

if L=[] then
 RETURN([]):
fi:

if nops(L)=1 then
 RETURN([1$L[1]]):
fi:

k:=nops(L):

for i from 1 to nops(L) while L[i]>1 do od:

L1:=[seq(L[i1]-1,i1=1..i-1)]:

[k,op(Conj(L1))]:
end:



#KickOut(w,i,j,n): Inputs a word w in {1,...,n-1}, kicks out w[i]...w[j] reananges it, and outputs the resulting permutation
#and its Sh1, and Sh2. try:
#KickOut([1, 2, 1, 4, 3, 2, 4, 1, 3, 2],1.3,5);
KickOut:=proc(w,i,j,n) local w1,pi1:
w1:=[op(j+1..nops(w),w),op(1..i-1,w)]:
pi1:=EvalMila(w1,n):
[pi1,[Sh1(pi1),Sh2(pi1)]]:
end:


#ReinerTable(w): Inputs a permutation w (of length n=nops(w)) and outputs the list 0f lists of length n-2 where each list has length inv(w)-2 , let's call it L such that
#L[j][k] is the set of minimal words that evaluate to w that have either [j,j+1,j] or [j+1,j,j+1] starting at the k-th place. Try:
#ReinerTable([5,4,3,2,1]);
ReinerTable:=proc(w) local n,gu,j,k,j1,k1,L,T,lu,gu1:
option remember:
 n:=nops(w):
 gu:=Milim(w):
 L:=inv(w):



for j from 1 to n-2 do
 for k from 1 to L-2 do
  T[j,k]:={}:
 od:
od:



for gu1 in gu do


 lu:=YBp(gu1):


 for k1 in lu[1] do
  j1:=gu1[k1]:
  T[j1,k1]:=T[j1,k1] union {gu1}:
 od:


 for k1 in lu[2] do
  j1:=gu1[k1+1]:
  T[j1,k1]:=T[j1,k1] union {gu1}:
 od:

od:



[seq([seq(T[j1,k1],k1=1..L-2)],j1=1..n-2)]:

end:






#CycS1(w,n): The cyclic shift of the word w to the left by one place, in the alphabet {1,...,n-1}. Try: 
#CycS1([1,2,1],3);
CycS1:=proc(w,n):[op(2..nops(w),w),n-w[1]]:end:


#CycS(w,n,r): The cyclic shift of the word w to the left by r places, in the alphabet {1,...,n-1}. Try: 
#CycS([1,2,1],3,2);
CycS:=proc(w,n,r) local gu,i:
gu:=w:
for i from 1 to r do
 gu:=CycS1(gu,n):
od:
gu:
end:




#Xnij(n,i,j): Inputs a positive integer n, and positive integer i between 1 and binomial(n,2)-2 and j such that 1<=j<=n-2, outputs The set of reduced words w of [n,...,1] that have the property 
#w[i]=j,w[i+1]=j+1, w[i+2]=j or w[i]=j+1,w[i+1]=j, w[i+2]=j+1 . Try:
#Xnij(5,1,1);
Xnij:=proc(n,i,j):

if not (1<=i and i<=binomial(n,2)-2 and 1<=j and j<=n-2) then
 RETURN({}):
fi:
ReinerTablePC(n)[j][i]:

end:



#Nus(n): The cardinalities that show up in the cardinalities of the intersections of Xnij(n1,i1,j1) and Xnij(n2,i2,j2). Try
#Nus(4)
Nus:=proc(n) local i1,j1,i2,j2,lu,gu:

gu:={}:

for i1 from 1 to binomial(n,2)-2 do
for i2 from i1+1 to binomial(n,2)-2 do
 for j1 from 1 to n-2 do
  for j2 from 1 to n-2 do
  lu:=nops(Xnij(n,i1,j1) intersect Xnij(n,i2,j2)):
  gu:=gu union {lu}:
 od:
od:
od:
od:
gu:
end:



#w0n(n): The permutation [n,n-1,...,1]
#w0n(10);
w0n:=proc(n) local i:
[seq(n+1-i,i=1..n)]:
end:


#w0nj(n,j): The permutation (j+2,j+1,j)*[n,n-1,...,1] in other words [n,n-1,..., j+3,j,j+1,j+2,j-1,...,1]. Try:
#w0nj(10,2);
w0nj:=proc(n,j) local i:
if not (type(n,integer) and type(j,integer) and n>=3 and j>=1 and j<=n-2) then
 print(`bad input`):
 RETURN(FAIL):
fi:
[seq(n+1-i,i=1..n-2-j),j,j+1,j+2,seq(n+1-i,i=n+2-j..n)]:
end:

#Shn(n): the partition [n-1,...,1]. Try:
#Shn(8);
Shn:=proc(n) local i:
[seq(n-i,i=1..n-1)]:
end:

#Shnj(n,j): the partition [n-1,...,1]- [0^(j-1),2,1,0^(n-2-j)]. It is the shape of w0nj(n,j).  Try:
#Shnj(10,1);
Shnj:=proc(n,j) local i:
if j<n-2 then
Shn(n) - [0$(j-1),2,1,0$(n-2-j)]:
else
[seq(n-i,i=1..n-3)]:
fi:

end:


#EXn(n): The expected number of Yang-Baxter moves on the symmetric group on n elements. Try:
#EXn(8);
EXn:=proc(n) local j,l,tot:
l:=binomial(n,2):
tot:=YF(Shn(n)):
(l-2)*2*add(YF(Shnj(n,j))/tot,j=1..n-2):
end:


#Xnrj1j2(n,r,j1,j2): The set of minimal words of w0n(n) with Yang-Baxter moves at location 1 of type j1, and at location r of type j2. Try:
#Xnrj1(6,3,1,2);
Xnrj1j2:=proc(n,r,j1,j2)

if not (r<=binomial(n,2)-2 and 1<=j1 and j1<=n-2 and 1<=j2 and j2<=n-2) then
 RETURN(FAIL):
fi:
Xnij(n,1,j1) intersect Xnij(n,r,j2):

end:


#Pars1(n,k): the set of partitions of n with largest part k. Try:
#Pars1(5,2);
Pars1:=proc(n,k) local k1,gu,gu1,gu11:
option remember:

if k>n then
 RETURN({}):
elif k=n then
 RETURN({[k]}):
else
gu:={}:

for k1 from 1 to min(k,n-k) do

 gu1:=Pars1(n-k,k1):
 gu:=gu union {seq([k,op(gu11)],gu11 in gu1)}:
od:
RETURN(gu):
fi:

end:


#Pars(n): all the partitions of n
Pars:=proc(n) local k:
{seq(op(Pars1(n,k)),k=1..n)}:
end:

#ParsL(L): The set of sub-shapes of L. Try: ParsL([5,4,3,2,1]);
ParsL:=proc(L) local gu,n,n1,gu1,gu11,gu11a,i1:

n:=convert(L,`+`):

gu:={}:
for n1 from 1 to n-1 do
gu1:=Pars(n1):
for gu11 in gu1 do
if nops(gu11)<=nops(L) then
gu11a:=[op(gu11),0$(nops(L)-nops(gu11))]:
if min(seq(L[i1]-gu11a[i1],i1=1..nops(L)))>=0 then
 gu:=gu union {gu11}:
fi:
fi:
od:
od:
gu:
end:

#Locatek(T,k): the places where k lives
Locatek:=proc(T,k) local i,j:
for i from 1 to nops(T) do
for j from 1 to nops(T[i]) do
 if T[i][j]=k then
  RETURN([i,j]):
 fi:
od:
od:
FAIL:
end:


#AllYF(n): All the YF(L) for L in the stair-case shape [n-1,n-2,...,1]. Try:
#AllYF(4);
AllYF:=proc(n) local gu,i,gu1:
gu:=ParsL([seq(n-i,i=1..n-1)]):
{seq(YF(gu1),gu1 in gu)}:
end:

#EGr(T,n): The reverse of the Edelman-Greene mapping
EGr:=proc(T,n) local w,i,j,k,lu,pi1,i1:

pi1:=[seq(i,i=1..n)]:

if nops(T)<>n-1 and {seq(seq(T[i][j],j=1..n-i),i=1..n-1)}<>{seq(i,i=1..n*(n-1)/2)} then
 RETURN(FAIL):
fi:

w:=[]:

for k from 1 to n*(n-1)/2 do


lu:=Locatek(T,k):


i:=n+1-lu[1]:
j:=lu[2]:

for i1 from 1 to n while pi1[i1]<>i do od:
w:=[op(w),i1-1]:
if pi1[i1-1]<>j then
 RETURN(FAIL):
fi:

pi1:=[op(1..i1-2,pi1),i,j,op(i1+1..n,pi1)]:
od:
w:
end:


#jeuT(T,n,k): Applying jeu-de-taquin to a standard tableu T of k+1,...,k+n . Try:
#jeuT([[1,3,4],[2,6],[5]]);

jeuT:=proc(T,n,k) local lu,T1,kha,i,j:
lu:=Locatek(T,n+k):

T1:=T:

while lu<>[1,1] do

i:=lu[1]:
j:=lu[2]:

 if i>1 and j>1 then

 kha:=max(T1[i-1][j], T1[i][j-1]):
 T1:=[op(1..i-1,T1), [op(1..j-1,T1[i]),kha,op(j+1..nops(T1[i]),T1[i])],op(i+1..nops(T1),T1)]:

  if T1[i-1][j]>T1[i][j-1] then
  lu:=[i-1,j]:
   else
  lu:=[i,j-1]:
 fi:

elif i=1 and j>1 then
 kha:= T1[i][j-1]:
 T1:=[[op(1..j-1,T1[1]),kha,op(j+1..nops(T1[1]),T1[1])],op(2..nops(T1),T1)]:
 lu:=[1,j-1]:

 elif i>1 and j=1 then
 kha:= T1[i-1][j]:
 T1:=[op(1..i-1,T1), [kha,op(2..nops(T1[i]),T1[i])],op(i+1..nops(T1),T1)]:
 lu:=[i-1,1]:
fi:
od:

[[k,op(2..nops(T1[1]),T1[1])],op(2..nops(T1),T1)]:

end:






#SchP(T,n): applies the Schutzenberger transform to tableu T of {1,...,n}. Try:
#SchP([[1,3,4],[2,6],[5]]);
SchP:=proc(T,n) local i,j,k,T1:

T1:=T:
for k from 0 to n-1 do
T1:=jeuT(T1,n,-k):
od:

[seq([seq(T1[i][j]+n,j=1..nops(T1[i]))],i=1..nops(T1))]:
end:



#TtoW(T,n): inputs a Young tableau of staircase shape [n-1,...,1] outputs a word in {1,..,n-1} that is a minimal decomposition of the permutation [n,n-1,..,1]. Try:
#TtoW([[1,2,3,4],[5,6,7],[8,9],[10]],5);
TtoW:=proc(T,n) local N,T1,k,w:
N:=n*(n-1)/2:
w:=[]:
T1:=T:


for k from 0 to N-1 do
 w:=[op(w),Locatek(T1,N-k)[2]]:
 T1:=jeuT(T1,N,-k):

od:

w:
end:





#MilimS(n): all the words representing minimal factorizations of [n,n-1,..,2,1], using TtoW(T,n). Try:
#MilimS(5);
MilimS:=proc(n) local gu,i,gu1:

gu:=SYT([seq(n-i,i=1..n-1)]):

{seq(TtoW(gu1,n),gu1 in gu)}:
end:


#RandW(n): a uniformly at random word representing [n,n-1,...,1]. Try:
#RandW(10);
RandW:=proc(n)  local i:
TtoW(GNW([seq(n-i,i=1..n-1)]),n):
end:



#Simu(n,K,L): generates K random maximal words  representing [n,n-1,...,1] and retruns the list whose i-th item is the ratio of those that have i-1 Yang-Baxter moves (up to the max), followed
#by a list of length L whose first item is the average, second the variance, followed by the moments about the mean. Try:
#Simu(10,100,4);
Simu:=proc(n,K,L) local i,x,f,M1,M2:

f:=evalf(add(x^YBtot(RandW(n)),i=1..K)/K):
M1:=[seq(coeff(f,x,i),i=0..degree(f,x))]:
M2:=evalf(AveAndMoms(f,x,L)):
M1,M2:

end:

#GNWs(n): A random Young tableau of staircase shape [n-1,...,1] using the Greene-Nijenhuis-Wilf algorithm. Try:
#GWNs(10);
GNWs:=proc(n) local i:GNW([seq(n-i,i=1..n-1)]):end:

#Ins11(row1,k): Given a list of increasing integers, row1, inserts
#k in the proper place of row1, and returns the bumped
#element, or 0, if it is at the end
Ins11:=proc(row1,k)
local i,m:
m:=nops(row1):
for i from 1 to m while row1[i]<=k do
od:

if i=m+1 then
RETURN([op(row1),k],0):
else

if not (row1[i]-k=1 and member(k,{op(1..i-1,row1)})) then
RETURN([op(1..i-1,row1),k,op(i+1..m,row1)],row1[i]):
else
RETURN([op(1..i-1,row1),k+1,op(i+1..m,row1)],row1[i]):
fi:

fi:

end:


#Ins1(tab1,k): Given a tableau tab1, inserts the element k
#and returns the modified tableau and the row number
#where it landed
Ins1:=proc(tab1,k)
local j,lu,k1,newrow,tab2:

k1:=k:
tab2:=tab1:
for j from 1 to nops(tab1) while k1<>0 do
lu:=Ins11(tab1[j],k1):
k1:=lu[2]:
newrow:=lu[1]:
tab2:=[op(1..j-1,tab2),newrow,op(j+1..nops(tab2),tab2)]:
k1:=lu[2]:
od:

if k1<>0 then
tab2:=[op(tab2),[k1]]:
RETURN(tab2,nops(tab2)):
else
RETURN(tab2,j-1):
fi:

end:


RS:=proc(perm)
local tab1,tab2,p,i,n,gu:
n:=nops(perm):
tab1:=[]: tab2:=[]:
for i from 1 to n do
gu:=Ins1(tab1,perm[i]):
tab1:=gu[1]:
p:=gu[2]:
if p<=nops(tab2) then
tab2:=[op(1..p-1,tab2),[op(tab2[p]),i],op(p+1..nops(tab2),tab2)]:
else
tab2:=[op(tab2),[i]]:
fi:

od:
tab1,tab2:
end:



#WtoT(w,n): inputs a tableau T of staircase  shape [n-1,...,1] and outputs the corresponding word representing a minimal factorization of the permutation [n,n-1,...,1]. According to
#Edelman-Greene. Try:
#WtoT([1,2,1],3);
WtoT:=proc(w,n) local i,N:
N:=n*(n-1)/2:
if nops(w)<>N then
 RETURN(FAIL):
fi:

RS([seq(w[N+1-i],i=1..N)])[2]:

end:


#PrXnij(n,i,j): Prints out the Tableaux images under the Edelman-Greene mapping of the words given by Xnij(n,i,j) (q.v.). Try
#PrXnij(4,1,1);
PrXnij:=proc(n,i,j) local gu,i1:
gu:=Xnij(n,i,j):

for i1 from 1 to nops(gu) do
PrT(WtoT(gu[i1],n)):
print(``):
od:


end:





#PrXnrj1j2(n,i,j): Prints out the Tableaux images under the Edelman-Greene mapping of the words given by Xnrj1j2(n,r,j1,j2) (q.v.). Try
#PrXnrj1j2(5,3,1,2);
PrXnrj1j2:=proc(n,r,j1,j2) local gu,i1:

gu:=Xnrj1j2(n,r,j1,j2):

for i1 from 1 to nops(gu) do
PrT(WtoT(gu[i1],n)):
print(``):
od:


end:



#DecXnijSlow(n,i,j): The set of pairs [pi1,pi2] such that inv(pi1)=i-1, inv(pi2)=n*(n-1)/2-i-2 and [n,...,1]=pi1 [j+2,j+1,j] pi2, via Xnij(n,i,j). Try:
#DecXnijSlow(5,1,1);
DecXnijSlow:=proc(n,i,j) local gu,gu1:
if not (type(n,integer) and n>2 and j>=1 and j<=n-2 and i>=1 and i<=n*(n-1)/2-2) then
print(`bad input `):
 RETURN(FAIL):
fi:
gu:=Xnij(n,i,j):
{seq([EvalMila([op(1..i-1,gu1)],n), EvalMila([op(i+3..nops(gu1),gu1)],n)],gu1 in gu)}:
end:

#PermuteList(n): The list of length binomial(n,2)+1 such that for i=0..binomial(n,2), L[i+1] is the set of permutations of length n with i inversions. Try:
#PermuteList(5);
PermuteList:=proc(n) local gu,i,gu1,T:
option remember:
gu:=convert(permute(n),set):

for i from 0 to binomial(n,2) do
T[i]:={}:
od:

for gu1 in gu do
 T[inv(gu1)]:=T[inv(gu1)] union {gu1}:
od:

[seq(T[i],i=0..binomial(n,2))]:
end:





#DecXnij(n,i,j): The set of pairs [pi1,pi2] such that inv(pi1)=i-1, inv(pi2)=n*(n-1)/2-i-2 and [n,...,1]=pi1 [j+2,j+1,j] pi2, via Xnij(n,i,j). The fast way. Try:
#DecXnij(5,1,1);
DecXnij:=proc(n,i,j) local GU1,GU2,Pij,W0,gu1,gu2,LU,i1,j1:
if not (type(n,integer) and n>2 and j>=1 and j<=n-2 and i>=1 and i<=n*(n-1)/2-2) then
 RETURN({}):
fi:

GU1:=PermuteList(n)[i]:
GU2:=PermuteList(n)[binomial(n,2)-i-1]:


Pij:=[seq(j1,j1=1..j-1),j+2,j+1,j,seq(j1,j1=j+3..n)]:
W0:=[seq(n+1-i1,i1=1..n)]:

LU:={}:

for gu1 in GU1 do
for gu2 in GU2 do
 if Mul(Mul(gu1,Pij),gu2)=W0 then
  LU:=LU union {[gu1,gu2]}:
 fi:
od:
od:

LU:

end:


#NuXnij(n,i,j): Same as nops(Xnij(n,i,j)) but done cleverly. Try:
#NuXnij(5,1,4);
NuXnij:=proc(n,i,j) local gu,gu1:
gu:=DecXnij(n,i,j):

2*add(nops(Milim(gu1[1]))*nops(Milim(gu1[2])),gu1 in gu):
end:

#EXnD(n): Like EXn(n) but done directly. Only for checiking.
EXnD:=proc(n) local j:
(binomial(n,2)-2)*add(NuXnij(n,trunc(binomial(n,2)/2),j),j=1..n-2)/YF([seq(n-j,j=1..n-1)]):
end:





#Xni1j1i2j2(n,i1,j1,i2,j2): The set of minimal words of w0n(n) with Yang-Baxter moves at location i1 of type j1, and at location i2 of type j2. Try:
#Xni1j1i2j2(6,1,1,9,1);
Xni1j1i2j2:=proc(n,i1,j1,i2,j2) local gu1,gu2:

Xnij(n,i1,j1) intersect Xnij(n,i2,j2):

end:


#LocatPlaces(T,S): Given a tableau T and a set of entries S, outputs the set of places where they reside. Try:
#LocatePlaces([[1,2,3,4],[5,6,7],[8,9]],{1,6});
LocatePlaces:=proc(T,S) local gu,i,j:
gu:={}:

for i from 1 to nops(T) do
 for j from 1 to nops(T[i]) do
  if member(T[i][j],S) then
   gu:=gu union {[i,j]}:
  fi:
 od:
od:
gu:
end:




#NuMilim(pi): the cardinality of Milim(pi) (q.v)
NuMilim:=proc(pi) local n,i,pi1,mu:
option remember:

if inv(pi)=0 then
 RETURN(1):
fi:

n:=nops(pi):

mu:=0:

for  i from 1 to n-1 do
  if pi[i]>pi[i+1] then
 pi1:=[op(1..i-1,pi),pi[i+1],pi[i],op(i+2..n,pi)]:
    mu:=mu+NuMilim(pi1):
  fi:
od:
mu:


end:

#InvTab(pi): The inversion Table  of the permutation pi. Try:
#InvTab([5,1,2,3,4]);
InvTab:=proc(pi) local n,co,i,j,gu:
n:=nops(pi):
gu:=[]:


for i from 1 to n do
co:=0:
 for j from i+1 to n do
 if pi[i]>pi[j] then
   co:=co+1:
 fi:
 od:
gu:=[op(gu),co]:
od:
gu:
end:



#MilaToCh(W,n): converts the word W in S_n to a chain. Try:
#MilaToCh([1,2,1],3);
MilaToCh:=proc(W,n) local i:
[seq(EvalMila([op(i..nops(W),W)],n),i=1..nops(W)+1)]:
end:


#tin(i,n): the transposition (i,i+1) of {1,...,n}. Try: tin(6,2);
tin:=proc(i,n) local j: 
if not (type(n,integer) and type(i,integer) and n>0 and i>0 and i<n) then
RETURN(FAIL):
fi:

[seq(j,j=1..i-1),i+1,i,seq(j,j=i+2..n)]:
end: