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



print(`First Written: June 2022: tested for Maple 2020 `):
print(`Version  : June 2022  `):
print():
print(`This is VicR.txt, one of the  Maple packages that`):
print(`accompany 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/YBM.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():





ezraC:=proc()
if args=NULL then

print(`The Checking procedures are: CheckRDgf, CheckRDnuij, CheckRDsetij, RDgfD, RDsetijD `):
print(``):

else
ezra(args):

fi:


end:

ezra1:=proc()
if args=NULL then

print(`The SUPPORTING procedures are: Des,  InvTab, IsVex, RDgfij, RDnuij, RDsetij `):
print(``):

else
ezra(args):

fi:


end:

ezra:=proc()
if args=NULL then

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

print(`The MAIN procedures are:  GFstory, RDgf,  RDnu, RDnx, RDset, YF `):
print(``):

elif nargs=1 and args[1]=CheckRDgf then
print(`CheckRDgf(n): checks that RDgf(pi,x)=RDgfD(pi,x) for all permutations of length n. Try:`):
print(`CheckRDgf(5);`):

elif nargs=1 and args[1]=CheckRDnuij then
print(`CheckRDnuij(w): checks RDnuij(w,i,j). Try:`):
print(`CheckRDnuij([5,4,3,2,1]);`):


elif nargs=1 and args[1]=CheckRDsetij then
print(`CheckRDsetij(n): Checks RDsetij(w,i,j) for all permutations w  of length n and all i,j in {1,..., n-1}. Try:`):
print(`CheckRDsetij(5);`):

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

elif nargs=1 and args[1]=GFstory then
print(`GFstory(N,x): The first N generating functions of weight-enumerators of the set of reduced words of [n,n-1,...,1] according to the statistic "number of Yang-Baxter moves" for n=1...N`):
print(`together with their averages and variances. Try:`):
print(`GFstory(7,x);`):


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]=IsVex then
print(`IsVex(pi): Is the permutation pi Vexilary (i.e. does it avoid the pattern 2143?. Try: `):
print(`IsVex([2,1,4,3]);`):

elif nargs=1 and args[1]=RDgf then
print(`RDgf(w,x): The generating function, in the avariable x,  of the set of reduced decomposiitons of the permutation w according to the statistic "number of Yang-Baxter moves". It uses dynamical programming. Try:`):
print(`RDgf([4,3,2,1],x);`):


elif nargs=1 and args[1]=RDgfD then
print(`RDgfD(w,x): The generating function, in the avariable x,  of the set of reduced decomposiitons of the permutation w according to the statistic "number of Yang-Baxter moves". Done directly. ONLY FOR CHECKING. Try:`):
print(`RDgfD([4,3,2,1],x);`):

elif nargs=1 and args[1]=RDgfij then
print(`RDgfij(w,i,j,x): The generating function according to the statistic "number of Yang-Baxter moves"  of reduced decomposiitons of the permutation w that start with [i,j]. Try:`):
print(`RDgfij([5,4,3,2,1],1,2,x);`):


elif nargs=1 and args[1]=RDnuij then
print(`RDnuij(w,i,j): The number of reduced decomposiitons of the permutation w that start with [i,j]. Try:`):
print(`RDnuij([5,4,3,2,1],1,2);`):

elif nargs=1 and args[1]=RDnu then
print(`RDnu(w): The number of reduced decomposiitons of the permutation w, using dynamical programming. Try:`):
print(`RDnu([4,3,2,1]);`):

elif nargs=1 and args[1]=RDnx then
print(`RDnx(n,x): The generating function, according to the statistic "number of Yang-Baxter moves" of the set of reduced permutations of [n,n-1,n-2,...,1]. `):
print(`In other words, RDgf([seq(n+1-i,i=1..n)],x); Try:`):
print(`[seq(RDnx(n,x),n=1..6)];`):

elif nargs=1 and args[1]=RDset then
print(`RDset(w): The set of reduced decomposiitons of the permutation w, using dynamical programming. Try:`):
print(`RDset([4,3,2,1]);`):

elif nargs=1 and args[1]=RDsetij then
print(`RDsetij(w,i,j): The set of reduced decomposiitons of the permutation w that start with [i,j].Done by Dynamical programming. Try:`):
print(`RDsetij([5,4,3,2,1],1,2);`):


elif nargs=1 and args[1]=RDsetijD then
print(`RDsetijD(w,i,j): The set of reduced decomposiitons of the permutation w that start with [i,j]. Done directly, for checking purposes only.Try:`):
print(`RDsetijD([5,4,3,2,1],1,2);`):



elif nargs=1 and args[1]=YF then
print(`YF(L): The number of Standard Young tableaux of shape L, using the Yang-Frobenius formula. Try: `):
print(`YF([3,2,1]);`):

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

fi:


end:


with(combinat):


##### Start From RedDec.txt
#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:

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

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


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

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


##### End From RedDec.txt


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

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

S:


end:

#RDset(w): The set of reduced decomposiitons of the permutation w,using dynamical programming. Try:
#RDset([4,3,2,1]);
RDset:=proc(w) local n,gu,lu,i,w1,gu1,gu11:
option remember:

n:=nops(w):

lu:=Des(w):

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


gu:={}:

for i in lu do
 w1:=[op(1..i-1,w),w[i+1],w[i],op(i+2..n,w)]:
 gu1:=RDset(w1):
 gu:=gu union {seq([i,op(gu11)],gu11 in gu1)}:
od:

gu:

end:



#RDnum(w): The number of reduced decomposiitons of the permutation w
RDnu:=proc(w) local n,gu,lu,i,w1,gu1:
option remember:

n:=nops(w):

lu:=Des(w):

if lu={} then
 RETURN(1):
fi:


gu:=0:

for i in lu do
 w1:=[op(1..i-1,w),w[i+1],w[i],op(i+2..n,w)]:
 gu1:=RDnu(w1):
 gu:=gu+gu1:
od:

gu:

end:




#RDsetij(w,i,j): The set of reduced decomposiitons of the permutation w that start with [i,j]
RDsetij:=proc(w,i,j) local w1,k,w2,lu,n,gu,gu1,gu11:
option remember:

if i=j then 
 RETURN({}):
fi:

n:=nops(w):

if inv(w)<2 then
 RETURN({}):
fi:

if inv(w)=2 then
 if EvalMila([i,j],n)=w then
  RETURN({[i,j]}):
 else
  RETURN({}):
 fi:
fi:



w2:=Mul(EvalMila([j,i],n),w):

if inv(w2)<>inv(w)-2 then
 RETURN({}):
fi:

w1:=Mul(EvalMila([i],n),w):

lu:=Des(w2):

gu:={}:

for k in lu do
 gu1:=RDsetij(w1,j,k):
 for gu11 in gu1 do
     gu:=gu union {[i,op(gu11)]}:
 od:
od:

gu:

end:


#RDnuij(w,i,j): The set of reduced decomposiitons of the permutation w that start with [i,j]
RDnuij:=proc(w,i,j) local w1,k,w2,lu,n,gu:
option remember:

if i=j then 
 RETURN(0):
fi:

n:=nops(w):

if inv(w)<2 then
 RETURN(0):
fi:

if inv(w)=2 then
 if EvalMila([i,j],n)=w then
  RETURN(1):
 else
  RETURN(0):
 fi:
fi:



w2:=Mul(EvalMila([j,i],n),w):

if inv(w2)<>inv(w)-2 then
 RETURN(0):
fi:

w1:=Mul(EvalMila([i],n),w):

lu:=Des(w2):

gu:=0:

for k in lu do
gu:=gu+RDnuij(w1,j,k):
od:

gu:

end:





#RDgfij(w,i,j,x): The set of reduced decomposiitons of the permutation w that start with [i,j]
RDgfij:=proc(w,i,j,x) local w1,k,w2,lu,n,gu:
option remember:

if i=j then 
 RETURN(0):
fi:

n:=nops(w):

if inv(w)<2 then
 RETURN(0):
fi:

if inv(w)=2 then
 if EvalMila([i,j],n)=w then
  RETURN(1):
 else
  RETURN(0):
 fi:
fi:



w2:=Mul(EvalMila([j,i],n),w):

if inv(w2)<>inv(w)-2 then
 RETURN(0):
fi:

w1:=Mul(EvalMila([i],n),w):

lu:=Des(w2):

gu:=0:

for k in lu do
 if  i=k and abs(j-i)=1 then
  gu:=expand(gu+x*RDgfij(w1,j,k,x)):
  else
  gu:=gu+RDgfij(w1,j,k,x):
 fi:
od:

gu:

end:

#RDgf(w,x): The weight-enumerator, according to the number of Yang-Baxter moves, of the reduced decompositions of the permutation w
#Try:
#RDgf([4,3,2,1],x):
RDgf:=proc(w,x) local n,i,j:
n:=nops(w):
if inv(w)<2 then
 1:
else
add(add(RDgfij(w,i,j,x),i=1..n-1),j=1..n-1):
fi:
end:

#RDnx(n,x):RDgf([n,n-1,...,1],x);
RDnx:=proc(n,x) local i:
RDgf([seq(n+1-i,i=1..n)],x) :
end:



#CheckRDnuij(w): checks RDnuij(w,i,j). Try:
#CheckRDnuij([5,4,3,2,1]);
CheckRDnuij:=proc(w) local n, i,j:
n:=nops(w):
[add(add(RDnuij(w,i,j),i=1..n-1),j=1..n-1), RDnu(w)]:
end:




#RDsetijD(w,i,j): The set of reduced decomposiitons of the permutation w that start with [i,j], done directly, for checking only
RDsetijD:=proc(w,i,j) local  lu,lu1,gu:
option remember:

lu:=RDset(w):


gu:={}:

for lu1 in lu do
 if nops(lu1)>=2 and lu1[1]=i and lu1[2]=j then
   gu:=gu union {lu1}:
 fi:
od:

gu:


end:



#CheckRDsetij1(w)
CheckRDsetij1:=proc(w) local n, i,j:
n:=nops(w):
evalb({seq(seq(evalb(RDsetij(w,i,j)=RDsetijD(w,i,j)),i=1..n-1),j=1..n-1)}={true}):
end:



#CheckRDsetij(n): Checks RDsetij(w,i,j) for all permutations w  of length n and all i,j in {1,..., n-1}. Try:
#CheckRDsetij(5);
CheckRDsetij:=proc(n) local gu,w:
gu:=permute(n):
evalb({seq(CheckRDsetij1(w),w in gu)}={true}):
end:


#IsVex(pi): Is the permutation pi Vexilary (i.e. does it avoid the pattern 2143?. Try
#IsVex([2,1,4,3]);
IsVex:=proc(pi) local n,i1,i2,i3,i4:

n:=nops(pi):

for  i1 from 1 to n-3 do
 for  i2 from i1+1 to n-2 do
  for  i3 from i2+1 to n-1 do
    for  i4 from i3+1 to n do
  if pi[i2]<pi[i1] and pi[i1]<pi[i4] and pi[i4]<pi[i3] then
    RETURN(false):
  fi:
  od:
 od:
od:
od:
true:
end:






#YB(w): The number of Yang-Baxter moves in a word w of {1,...,n-1}
YB:=proc(w) local co,i:

co:=0:

for i from 1 to nops(w)-2 do

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


co:

end:

#RDgfD(pi,x): The generating function of the reduced words that evaluate to pi according to the weight x^YB(w) done DIRECTLY. Try:
#RDgfD([4,3,1,2],x);
RDgfD:=proc(pi,x) local gu,w:

gu:=RDset(pi):


add(x^YB(w),w in gu):

end:

#CheckRDgf(n): checks that RDgf(pi,x)=RDgfD(pi,x) for all permutations of length n. Try:
#CheckRDgf(5);
CheckRDgf:=proc(n) local gu,pi,x:
gu:=permute(n):
evalb({seq(RDgf(pi,x)-RDgfD(pi,x),pi in gu)}={0}):
end:

#GFstory(N,x): The first N generating functions of weight-enumerators of the set of reduced words of [n,n-1,...,1] according to the statistic "number of Yang-Baxter moves" for n=1...N
#together with their averages and variances. Try:
#GFstory(7,x);
GFstory:=proc(N,x) local n,gu,gu1,i:

gu:=[]:

for n from 1 to N do
gu1:=RDnx(n,x):
print(`The generating function of the number of reduced decompostions of `, [seq(n+1-i,i=1..n)], ` is `, gu1):
print(``):
print(`and in Maple notation `):
print(``):
lprint(gu1):
print(``):
print(`The expectation is`, subs(x=1,diff(gu1,x))/subs(x=1,gu1), `and the variance is`,  subs(x=1,diff(gu1,x,x))/subs(x=1,gu1)):
print(``):
gu:=[op(gu),gu1]:
print(``):
od:

print(`To sum up, the first `,N, `generating functions are `):
print(``):
lprint(gu):
print(``):
gu:
end: