######################################################################
##MahonianStat: Save this file as MahonianStat                       #
#To use it, stay in the                                              #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read MahonianStat : <Enter>                        #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Rutgers University ,                  #
#zeilberg at math dot rutgers dot edu                                #
######################################################################
 
#Created: Aug. 2008
 
print(`Created: Aug. 2008`):
print(` This is MahonianStat `):
print(`to  prove asymptotic normality and compute`):
print(`symbolically(!) moments of the r.v. "number of inversions" `):
print(`in the set of words with a[1] 1's, ..., a[n] n's `):
print(`in terms of explicit polynomials in the symmetric`):
print(`functions of the a[1], ..., a[n] `):
print(`for any desired (numeric) moments, `):
print(`and, for the two-letter case, even for symbolic moment,`):
print(` asymptotically, to any desired order.`):
print(``):
print(`written by Doron Zeilberger.`):
print(``):
print(`It  accompanies the article `):
print(`The Mahonian Probability Distribution`):
print(`on Words is Asymptotically Normal`):
print(`by Rodney Canfield, Svante Janson and Doron Zeilberger`):
print(`available from Zeilberger's website (and possibly other places)`):
print(``):
print(`Please report bugs to zeilberg at math dot rutgers dot edu`):
print(``):
 print(`The most current version of this  package and paper`):
 print(` are  available from`):
 print(`http://www.math.rutgers.edu/~zeilberg/papers1.html  .`):
 print(`For a list of the procedures type ezra();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):

with(combinat):

ezraPerm1:=proc()

if args=NULL then
 print(` The supporting procedures are: CheckAppxInv`):
 print( ` EXijk, Fni, Fn, BMomAsyE, BMomAsyO,`):
 print(` GuessPolG, MXijk, Pele , StirlingS2 `):
else
ezra(args):
fi:

end:


ezraPerm:=proc()

if args=NULL then
 print(`The main procedures are: `):
 print(` Alpha, AlphaInv, AlphaInvAsy, AlphaInvAsyS`):
 print(`AppxInv, ithMomInv, ithBMomInv, MomAsy,`):
 print(`MXijkS `):
 print(` `):



elif nops([args])=1 and op(1,[args])=Alpha then
print(`Alpha(n,i,j,k): The generalized covariance of`):
print(`des,inv, maj, with powers i,j,k, divided by`):
print(`the (variances of des)^(i/2)*(variances of inv)^(i/2)*`):
print(`(variances of maj)^(i/2)`):
print(`For example, to get the Kurtosis of des`):
print(`try: Alpha(n,4,0,0);`):

elif nops([args])=1 and op(1,[args])=AlphaInv then
print(`AlphaInv(n,i): The alpha_i for the r.v. "no. of inversions"`):
print(`on S_n, for sybolic n and numeric i.`):
print(`For example, for the Kurtosis do: AlphaInv(n,4);`):


elif nops([args])=1 and op(1,[args])=AlphaInvAsyS then
print(`AlphaInvAsyS(n,r,s): Given symbolic n and r, and`):
print(`a pos. integer s, finds`):
print(`the first s terms in the asymptotic`):
print(`expansion of the parameter alpha_2r `):
print(`(the 2r-th moment about the mean divided `):
print(`by the r-th power of the variance`):
print(`for the r.v. number of inversions`):
print(`on n-perms. For example, try:`):
print(`AlphaInvAsyS(n,r,4);`):

elif nops([args])=1 and op(1,[args])=AlphaInvAsy then
print(`AlphaInvAsy(n,i,R): The asymptotics of`):
print(`alpha_i for the r.v. "no. of inversions"`):
print(`on S_n, for sybolic n and numeric i, to R terms.`):
print(`For example, for the asymptotics Kurtosis to 6 terms `):
print(`do: AlphaInvAsy(n,4,6);`):

elif nops([args])=1 and op(1,[args])=AppxInv then
print(`AppxInv(n,x): An approximtion for`):
print(`the number of n-perms with x inversions.`):
print(`For example, try: AppxInv(50,binomial(50,2)/2);`):



elif nops([args])=1 and op(1,[args])=BMomAsyE then
print(`BMomAsyE(n,r,s): inputs symbols n and r, and a`):
print(`a pos. integer s, and outputs a`):
print(`guess for the leading terms, of ithBMomInv(n,2*r) `):
print(`from n^(3*r), to n^(3*r-s)  in the expression`):
print(`for the 2r-th binomial moment of the r.v. "numer of inversions"`):
print(`acting on the set of n-perms. For example, try:`):
print(`BMomAsyE(n,r,4);`):

elif nops([args])=1 and op(1,[args])=BMomAsyO then
print(`BMomAsyO(n,r,s): inputs symbols n and r, and a`):
print(`a pos. integer s, and outputs a`):
print(`guess for the leading terms, of ithBMomInv(n,2*r+1) `):
print(`from n^(3*r), to n^(3*r-s)  in the expression`):
print(`for the (2r+1)-th binomial moment of the r.v. "numer of inversions"`):
print(`acting on the set of n-perms. For example, try:`):
print(`BMomAsyO(n,r,4);`):

elif nops([args])=1 and op(1,[args])=CheckAppxInv then
print(`CheckAppxInv(n1): cheks AppxInv(n,x)`):
print(`the number of n-perms with x inversions.`):
print(`For example, try: CheckAppxInv(50);`)

elif nops([args])=1 and op(1,[args])=CheckBetterAppxInv then
print(`CheckBetterAppxInv(n1,s): cheks BetterAppxInv(n,x,s)`):
print(`the number of n-perms with x inversions.`):
print(`For example, try: CheckBetterAppxInv(50,2);`):

elif nops([args])=1 and op(1,[args])=EXijk then
print(`EXijk(n,i,j,k): the sum of des(pi)^i*inv(pi)^j*maj(pi)^k`):
print(`over all n-perms divided by n!. For example, try`):
print(`EXijk(3,1,1,2);`):

elif nops([args])=1 and op(1,[args])=Fn then
print(`Fn(n,t,q1,q2): the weight enumerator, according to the weight`):
print(`t^des(pi)*q1^inv(pi)*q2^maj(pi) of all n-perms. `):
print(`For example, try`):
print(`Fn(3,t,q1,q2);`):

elif nops([args])=1 and op(1,[args])=Fni then
print(`Fni(n,i,t,q1,q2): the weight enumerator, according to the weight`):
print(`t^des(pi)*q1^inv(pi)*q2^maj(pi) of all n-perms that end in i.`):
print(`For example, try`):
print(`Fni(3,2,t,q1,q2);`):

elif nops([args])=1 and op(1,[args])=GuessPolG then
print(`GuessPolG(L,n): guesses a polynomial of degree d in n for`):
print(` the list L, such that P(i)=L[i] for i=s0..nops(L) for `):
print(`some s0. For example, try: `):
print(`GuessPolG([seq(i,i=1..10)],n);`):

elif nops([args])=1 and op(1,[args])=ithBMomInv then
print(`ithBMomInv(i,n): a polynomial expression in n, for `):
print(` E[binomial(X-av,r)] `):
print(` for the r.v. "number of inversions of n-perms `):
print(`where n is symbolic.`):
print(`For example, try:`):
print(`ithBMomInv(4,n);`):

elif nops([args])=1 and op(1,[args])=ithMomInv then
print(`ithMomInv(i,n): a polynomial expression in n, for the i-th moment`):
print(`about the mean for the number of inversions of n-perms `):
print(`where n is symbolic. The same as MXijkS(n,0,i,0), but faster.`):
print(`For example, try:`):
print(`ithMomInv(4,n);`):

elif nops([args])=1 and op(1,[args])=MomAsy then
print(`MomAsy(n,r,s): Given symbols n and r and a pos. integer s,`):
print(`the asympt. of the 2r-th moment`):
print(`fron n^(3*r) to n^(3*r-s)`):

elif nops([args])=1 and op(1,[args])=MXijk then
print(`MXijk(n,i,j,k): the sum of `):
print(`[des(pi)-av(des)]^i*[inv(pi)-av(inv)]^j*[maj(pi)-av(maj)]^k`):
print(`over all n-perms divided by n!. For example, try`):
print(`MXijk(3,1,1,2);`):

elif nops([args])=1 and op(1,[args])=MXijkS then
print(`MXijkS(n,i,j,k): conjectures a polynomial expression`):
print(`in the symbol n, for the sum of `):
print(`[des(pi)-av(des)]^i*[inv(pi)-av(inv)]^j*[maj(pi)-av(maj)]^k`):
print(`over all n-perms divided by n!. For example, try`):
print(`MXijkS(n,1,1,2);`):


elif nops([args])=1 and op(1,[args])=Pele then
print(`Pele(z,n,R):The first R terms of the formal power series`):
print(`((1+z)^n-(1+z)^(-n))/((1+z)-(1+z)^(-1));`):
print(`For example, try:`):
print(`Pele(z,n,5);`):



elif nops([args])=1 and op(1,[args])=StirlingS2 then
print(`StirlingS2(n,r): The polynomial that expresses Stirling2(n,n-r).`):
print(`For example, try: StirlingS2(n,2);`):

else
print(`There is no ezra for`,args):
fi:
 
end:


##begin guessing polynomials
#GuessPol1(L,d,n): guesses a polynomial of degree d in n for
# the list L, such that P(i)=L[i+1] for i=1..nops(L)
#For example, try: 
#GuessPol1([(seq(i,i=1..10),1,n);
GuessPol1:=proc(L,d,n) local P,i,a,eq,var:
if d>nops(L)-2 then
 ERROR(`the list is too small`):
fi:

P:=add(a[i]*n^i,i=0..d):
var:={seq(a[i],i=0..d)}:
eq:={seq(subs(n=i,P)-L[i],i=1..nops(L))}:

var:=solve(eq,var):

if var=NULL then
 RETURN(FAIL):
fi:

subs(var,P):

end:



#GuessPol(L,n): guesses a polynomial of degree d in n for
# the list L, such that P(i)=L[i] for i=1..nops(L) for example, try: 
#GuessPol([(seq(i,i=1..10),n);
GuessPol:=proc(L,n) local d,gu:

for d from 0 to nops(L)-2 do
 gu:=GuessPol1(L,d,n):
 if gu<>FAIL then
    RETURN(gu):
 fi:
od:

FAIL:

end:



#GuessPolG(L,n): guesses a polynomial of degree d in n for
# the list L, such that P(i)=L[i] for i=s0..nops(L) for 
#some s0. For example, try: 
#GuessPolG([seq(i,i=1..10)],n);
GuessPolG:=proc(L,n) local gu,s0,L1:

for s0 from 1 to nops(L)-4 do
L1:=[op(s0..nops(L),L)]:
gu:=GuessPol(L1,n):

 if gu<>FAIL then
  RETURN(expand(subs(n=n-s0+1,gu))):
fi:
od:
FAIL:
end:





#Fni(n,i,t,q1,q2): the weight enumerator, according to the weight
#t^des(pi)*q1^inv(pi)*q2^maj(pi) of all n-perms that end in i
#For example, try
#Fni(3,2,t,q1,q2);
Fni:=proc(n,i,t,q1,q2) local gu,j:
option remember:

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

gu:=0:

for j from 1 to i-1 do
 gu:=gu+Fni(n-1,j,t,q1,q2):
od:

for j from i to n-1 do
 gu:=gu+Fni(n-1,j,t,q1,q2)*t*q2^(n-1):
od:

expand(q1^(n-i)*gu):

end:





#Fn(n,t,q1,q2): the weight enumerator, according to the weight
#t^des(pi)*q1^inv(pi)*q2^maj(pi) of all n-perms 
#For example, try
#Fn(3,t,q1,q2);
Fn:=proc(n,t,q1,q2) local i: add(Fni(n,i,t,q1,q2),i=1..n):end:


#EXijk(n,i,j,k): the sum of des(pi)^i*inv(pi)^j*maj(pi)^k
#over all n-perms divided by n!. For example, try
#Exijk(3,1,1,2);
EXijk:=proc(n,i,j,k) local gu,t,q1,q2,i1:
option remember:
gu:=Fn(n,t,q1,q2):

for i1 from 1 to i do
 gu:=t*diff(gu,t):
od:

for i1 from 1 to j do
 gu:=q1*diff(gu,q1):
od:


for i1 from 1 to k do
 gu:=q2*diff(gu,q2):
od:

subs({q1=1,q2=1,t=1},gu)/n!:

end:



#MXijk(n,i,j,k): the sum of 
#[des(pi)-av(des)]^i*[inv(pi)-av(inv)]^j*[maj(pi)-av(maj)]^k
#over all n-perms divided by n!. For example, try
#Mxijk(3,1,1,2);
MXijkOld:=proc(n,i,j,k) local gu,a1,a2,a3,i1,j1,k1:
option remember:

a1:=EXijk(n,1,0,0):
a2:=EXijk(n,0,1,0):
a3:=EXijk(n,0,0,1):

gu:=0:

for i1 from 0 to i do
 for j1 from 0 to j do
  for k1 from 0 to k do
   gu:=gu+binomial(i,i1)*binomial(j,j1)*binomial(k,k1)*
   (-a1)^(i-i1)*(-a2)^(j-j1)*(-a3)^(k-k1)*EXijk(n,i1,j1,k1):
  od:
 od:
od:
gu:
end:




#MXijkS(n,i,j,k): Given a symbol n and non-neg. integers
#i,j,k,
#conjectures (a rigorous!) polynomial expression in n
#for the average of
#[des(pi)-av(des)]^i*[inv(pi)-av(inv)]^j*[maj(pi)-av(maj)]^k
#over all n-perms divided by n!. For example, try
#MxijkS(n,1,1,2);
MXijkS:=proc(n,i,j,k) local n0,gu:
gu:=[seq(MXijk(n0,i,j,k),n0=1..i+j+k+6)]:
factor(GuessPolG(gu,n)):
end:





###moments of Mahonian statistics
F:=proc(n,q) local i:
expand(normal(mul((1-q^i),i=1..n)/(1-q)^n)):
end:

#Fi(i,n): The i^th moment of the r.v. number of inversions
#over perms
Fi:=proc(i,n1) local j,gu,q:
gu:=F(n1,q)/q^(n1*(n1-1)/4):

for j from 1 to i do
 gu:=q*diff(gu,q):
od:
subs(q=1,gu)/n1!:
end:


#Gi(i,n): E[binomial(X,r)] of the r.v. number of inversions
#over perms
Gi:=proc(i,n1) local j,gu,q:
gu:=F(n1,q)/q^(n1*(n1-1)/4):

for j from 1 to i do
 gu:=diff(gu,q):
od:
subs(q=1,gu)/n1!/i!:
end:

#ithMomInvSlow(i,n): Like ithMomInv, but much slower, using
#polynomial fitting. For checking purposes only
ithMomInvSlow:=proc(i,n) local n1,N:
option remember:
N:=2*i+5:
factor(GuessPolG([seq(Fi(i,n1),n1=1..N)],n)):
end:

ithBMomInvSlow:=proc(i,n) local n1,N:
option remember:
N:=2*i+5:
factor(GuessPolG([seq(Gi(i,n1),n1=1..N)],n)):
end:


###end moments of Mahonian statistics

#MXijk(n,i,j,k): the sum of (des(pi)-av(des))^i*inv(pi)^j*maj(pi)^k
#over all n-perms divided by n!. For example, try
#Mxijk(3,1,1,2);
MXijk:=proc(n,i,j,k) local gu,t,q1,q2,i1:
option remember:
gu:=Fn(n,t,q1,q2)/t^((n-1)/2)/q1^(n*(n-1)/4)/q2^(n*(n-1)/4):

for i1 from 1 to i do
 gu:=t*diff(gu,t):
od:

for i1 from 1 to j do
 gu:=q1*diff(gu,q1):
od:


for i1 from 1 to k do
 gu:=q2*diff(gu,q2):
od:

subs({q1=1,q2=1,t=1},gu)/n!:

end:


Alpha:=proc(n,i,j,k)
MXijkS(n,i,j,k)/MXijkS(n,2,0,0)^(i/2)/MXijkS(n,0,2,0)^(j/2)
/MXijkS(n,0,0,2)^(k/2):
end:


#Pele(z,n,R):The first R terms of the formal power series
#((1+z)^n/2-(1+z)^(-n/2))/((1+z)^(1/2)-(1+z)^(-1/2));
#For example, try:
#Pele(z,n,5);
Pele:=proc(z,n,R) local gu,i:
gu:=((1+z)^(n/2)-(1+z)^(-n/2))/((1+z)^(1/2)-(1+z)^(-1/2))/n;
gu:=taylor(gu,z=0,R+3):
add(expand(coeff(gu,z,i))*z^i,i=0..R):
end:




#ithBMomInv(i,n): the polynomial expression,
#that describes  E[binomial(X-av,r)] where X is the
#r.v. "number of inversions" on the set of n-perms."
#and av is the average. For example, try:
#ithBMomInv(2,n);
ithBMomInv:=proc(r,n)
local gu,b,s,z,m:
option remember:
if r=0 then
 RETURN(1):
fi:

gu:=Pele(z,n,r):
b:=add(coeff(gu,z,s)*subs(n=n-1,ithBMomInv(r-s,n)),s=1..r);
b:=sum(b,n=1..m):
factor(subs(m=n,b)):
end:


#PowerToBinom(x,r): Given r, finds the vector of length
#r+1 such that x^r=a[0]+a[1]*binomial(x,1)+...+a[r]*binomial(x,r)
#For example, try:
#PowerToBinom(2);
PowerToBinom:=proc(r) local a,i,eq,var,gu,x:
var:={seq(a[i],i=0..r)}:
gu:=expand(x^r-expand(add(a[i]*binomial(x,i),i=0..r))):
eq:={seq(coeff(gu,x,i),i=0..r)}:
var:=solve(eq,var):
[seq(subs(var,a[i]),i=0..r)]:
end:


ithMomInv:=proc(i,n) local gu,r:
gu:=PowerToBinom(i):
factor(add(gu[r+1]*ithBMomInv(r,n),r=0..i)):
end:


#AlphaInv(n,i): The alpha_i for the r.v. "no. of inversions"
#on S_n, for sybolic n and numeric i.
#For example, for the Kurtosis do: AlphaInv(n,4);

AlphaInv:=proc(n,i)
normal(ithMomInv(i,n)/ithMomInv(2,n)^(i/2)):
end:


#AlphaInvAsy(n,i,R): The asymptotics of
#alpha_i for the r.v. "no. of inversions"
#on S_n, for sybolic n and numeric i, to R terms
#For example, for the asymptotics Kurtosis to 6 terms 
#do: AlphaInvAsy(n,4,6);
AlphaInvAsy:=proc(n,i,R) local gu,z,i1:
gu:=normal(ithMomInv(i,n)/ithMomInv(2,n)^(i/2)):
gu:=normal(subs(n=1/z,gu)):
gu:=taylor(gu,z=0,R+3):
add(coeff(gu,z,i1)/n^i1,i1=0..R):

end:

#BMomAsyE(n,r,s): inputs symbols n and r, and a
#a pos. integer K, and outputs a
#guess for the leading terms, of ithBMomInv(
#from n^(3*r), tp n^(3*r-s)  in the expression
#for the 2r-th binomial moment of the r.v. "numer of inversions"
#acting on the set of n-perms. For example, try:
#BMomAsyE(n,r,2);
BMomAsyE:=proc(n,r,s) local gu,i,mu,lu,lu1,K,i1:
K:=2*s+4:
mu:=[seq(ithBMomInv(2*i,n),i=1..K)]:

gu:=0:

for i from 0 do
 lu:=[seq(coeff(mu[i1],n,3*i1-i)*72^i1*i1!,i1=1..K)]:
 lu1:=GuessPol(lu,r):
 if lu1=FAIL then
   RETURN((1/72^r/r!)*gu):
 fi:
 gu:=gu+lu1*n^(3*r-i):
od:

end:




#BMomAsyO(n,r,s): inputs symbols n and r, and a
#a pos. integer s and outputs a
#guess for the leading terms, of ithBMomInv(n,2*r+1)
#from n^(3*r+1)  to n^(3*r+1-s)  in the expression
#for the 2r-th binomial moment of the r.v. "numer of inversions"
#acting on the set of n-perms. For example, try:
#BMomAsyE(n,r,2);
BMomAsyO:=proc(n,r,s) local gu,i,mu,lu,lu1,K,i1:
K:=2*s+4:
mu:=[seq(ithBMomInv(2*i+1,n),i=1..K)]:

gu:=0:

for i from 0 do
 lu:=[seq(coeff(mu[i1],n,3*i1-i)*72^i1*i1!,i1=1..K)]:
 lu1:=GuessPol(lu,r):
 if lu1=FAIL then
   RETURN((1/72^r/r!)*gu):
 fi:
 gu:=gu+lu1*n^(3*r-i):
od:

end:







#StirlingS2(n,r): The polynomial that expresses Stirling1(n,n-r)
#For example, try: StirlingS2(n,2);
StirlingS2:=proc(n,r) local gu,n1:

gu:=GuessPol([seq(stirling2(n1,n1-r),n1=1..2*r+6)],n):
if gu=FAIL then
 RETURUN(FAIL):
fi:
if {seq(subs(n=n1,gu)-stirling2(n1,n1-r),n1=1..4*r+6)}<>{0} then
  RETURN(FAIL):
else
 RETURN(gu):
fi:

end:



PowerToBinomS:=proc():
end:


#MomAsy(n,r,s): Given symbols n and r and a pos. integer s,
#the asympt. of the 2r-th moment
#fron n^(3*r) to n^(3*r-s)
MomAsy:=proc(n,r,s) local gu,i,s1,muE,muO,rosh:

muE:=BMomAsyE(n,r,s):
muO:=BMomAsyO(n,r,s):

gu:=
add(subs(r=2*r,StirlingS2(r,2*s1))*(2*r-2*s1)!*subs(r=r-s1,muE),s1=0..s)+
add(subs(r=2*r,StirlingS2(r,2*s1-1))*(2*r-2*s1+1)!*subs(r=r-s1,muO),s1=1..s):

rosh:=n^(3*r)*(2*r)!/72^r/r!:

gu:=expand(gu/rosh):

rosh*(add(factor(coeff(gu,n,-i))/n^i,i=0..s)):
end:





#AlphaInvAsyS(n,r,s): Given symbolic n and r, and
#a pos. integer s, finds
#the first s terms in the asymptotic
#expansion of the parameter alpha_2r 
#(the 2r-th moment about the mean divided 
#by the r-th power of the variance`):
#for the r.v. #of inversions
#on n-perms. For example, try:
#AlphaInvAsyS(n,r,4);
AlphaInvAsyS:=proc(n,r,s) local lu,i,x:

lu:=simplify(MomAsy(n,r,s)/n^(3*r)):
lu:=normal(lu/subs(r=1,lu)^r):
lu:=normal(subs(n=1/x,lu)):
lu:=normal(lu/((2*r)!/r!/2^r)):
lu:=taylor(lu,x=0,s+2):
lu:=add(eval(coeff(lu,x,i))/n^i,i=0..s):
lu:=(2*r)!/r!/2^r*(add(factor(simplify(coeff(lu,n,-i)))/n^i,i=0..s)):
end:




#AppxInv(n,x): An approximtion for
#the number of n-perms with x inversions.
#For example, try: AppxInv(50,binomial(50,2)/2);
AppxInv:=proc(n1,x) local y,sig,n:
sig:=evalf(sqrt(subs(n=n1,ithMomInv(2,n)))):
y:=evalf(binomial(n1,2)/2-x)/sig:
if abs(y)>2 then
 RETURN(FAIL):
fi:

evalf(n1!*exp(-y^2/2)/sqrt(2*Pi)/sig):

end:









CheckAppxInv:=proc(n1) local x,av,gu,q,sig,n,i:
av:=n1*(n1-1)/4:
sig:=evalf(sqrt(subs(n=n1,ithMomInv(2,n)))):
gu:=expand(mul(normal((1-q^i)/(1-q)),i=1..n1)):
{seq(evalf(AppxInv(n1,x)/coeff(gu,q,x)),x=trunc(av-sig)..trunc(av+sig))}:
end:




############End of permutation parts#####################

###########Begin words 


ezra1:=proc()

if args=NULL then
 print(` The supporting procedures are: AlphaW2t, AsyAlphaW2t, `):
 print(`,Bdok `):
 print(` CheckAppxW2, CheckBetterAppxW2, `):
 print(` CheckBetterAppxW2kFixed`):
 print(` CheckithBMk1,CompareAppxWk, Fabq,Gabq `):
 print(`ithBM2, ithBM3,ithBMk, ithMomW2t, NithMomW2t,`):
 print(`NithMomWkt, PeleW2, PeleWk  `):
 print(` PolToSym `):

else
ezra(args):
fi:

end:

ezra:=proc()

if args=NULL then
 print(`The main procedures are: AlphaWktE, AppxW2, AppxWk, AsyAlphaW2tS, `):
 print(` AsyAlphaWktE `):
 print(`  AsyFixedkAlpha, AsyMgfW2kFixed, AsyMgfW2t `):
 print(`AsyNithMomW2t, AsyPhiW2kFixed,`):
 print(`AsyPhiW2t, AsyRthMomW2t, BetterAppxW2, `):
 print(` BetterAppxW2kFixed, FixedkAlpha, ithMomW2, `):
 print(`ithMomWk, ithMomW2t , ithMomWktE, ithMomWktM`):


elif nops([args])=1 and op(1,[args])=AlphaW2t then
print(`AlphaW2t(r1,a,b,t): alpha_{2r1} of the r.v. Number of Inversions`):
print(`with a*t 1's and b*t 2's for numeric r1. For example, try:`):
print(`AlphaW2t(3,a,b,t);`):

elif nops([args])=1 and op(1,[args])=AlphaWktE then
print(`AlphaWktE(r1,e,t): alpha_{2r1} of the r.v. Number of Inversions`):
print(`with a[1]*t 1's, a[2]*t 2's ..., a[n]*t n's for numeric r1.`):
print(` For example, try:`):
print(`AlphaWktE(2,e,t);`):

elif nops([args])=1 and op(1,[args])=AppxW2 then
print(`AppxW2(a1,b1,x): An approximtion for`):
print(`the number of words in {1^a1,2^b1}  with x inversions.`):
print(`For example, try: AppxW2(20,20,200);`):

elif nops([args])=1 and op(1,[args])=AppxWk then
print(`AppxWk(L,x): An approximtion for`):
print(`the number of words in {1^L[1],2^L[2], ..., n^L[n]}`):
print(`  with x inversions. For example, try: `):
print(`AppxWk([20,20,20],600);`):

elif nops([args])=1 and op(1,[args])=AsyAlphaW2t then
print(`AsyAlphaW2t(r,a,b,t,s): The first s terms in the`):
print(`asymptotic expansion for alpha_{2r} (normlaized 2r-th moment)`):
print(`for the r.v. NumberOfInversions on words with a*t 1's and b*t 2's.`):
print(`for numeric r, but symbolic a,b,t`):
print(`For example, try:`):
print(`AsyAlphaW2t(4,a,b,t,2);`);

elif nops([args])=1 and op(1,[args])=AsyAlphaWktE then
print(`AsyAlphaWktE(r,e,t,s): The first s terms in the`):
print(`asymptotic expansion (in t) `):
print(`for alpha_{2r} (normlaized 2r-th moment)`):
print(`for the r.v. NumberOfInversions on words with1`):
print(` a[1]*t 1's a[2]*t 2's, ..., a[n]*t n's`):
print(`in terms the elementary symmetric functions e[1],e[2], ...`):
print(`of a[1], ..., a[n], and t`):
print(`for numeric r, but symbolic a,b,t`):
print(`For example, try:`):
print(`AsyAlphaWktE(4,e,t,2);`);

elif nops([args])=1 and op(1,[args])=AsyAlphaW2tS then
print(`AsyAlphaW2tS(r,a,b,t,s): The first s terms in the`):
print(`asymptotic expansion for alpha_{2r} (normlaized 2r-th moment)`):
print(`for the r.v. NumberOfInversions on words with a*t 1's and b*t 2's.`):
print(`Here r,a,b,t are symbols, but s is a numeric pos. integer`):
print(`For example, try:`):
print(`AsyAlphaW2tS(r,a,b,t,2);`);


elif nops([args])=1 and op(1,[args])=AsyFixedkAlpha then
print(`AsyFixedkAlpha(r,k,s): the asymptotic normalized 2i-th moment`):
print(`for the number of inversions with k 1's and n-k 2's`):
print(`with a fixed k, as n goes to infinity, up to s terms`):
print(`in 1/k. For example, try:`):
print(`AsyFixedkAlpha(r,k,2);`):


elif nops([args])=1 and op(1,[args])=AsyMgfW2kFixed then
print(`AsyMgfW2kFixed(k,x,s):`):
print(`The asymptotic moment generating function of the normalized`):
print(`r.v. "number of  inversions" of k 1's and n-k 2's up to order s`):
print(`with k fixed and n large`):
print(`in terms of x. For example, try:`):
print(`AsyMgfW2kFixed(k,x,2);`):

elif nops([args])=1 and op(1,[args])=AsyMgfW2t then
print(`The asymptotic moment generating function of the normalized`):
print(`r.v. "number of  inversions" of a*t 1's and b*t 2's up to order s`):
print(`in terms of x. For example, try: `):
print(`AsyMgfW2t(a,b,t,2,x);`):

elif nops([args])=1 and op(1,[args])=AsyNithMomW2t then
print(`AsyNithMomW2t(r,a,b,t,s): the first s terms in the`):
print(`asymptotic expansion, w.r.t. t, of the NORMALIZED 2r-th moment `):
print(`(about the mean) of the`):
print(`r.v. Number of inversions for words with a*t 1's and b*t 2's`):
print(`Here a,b,t, and r are symbolic(!), only s is numeric.`):
print(`For example, try:`):
print(`AsyNithMomW2t(r,a,b,t,3); `):




elif nops([args])=1 and op(1,[args])=AsyPhiW2kFixed then
print(`AsyPhiW2kFixed(k,2,x): the asympt. to order s of the`):
print(`continuous limit `):
print(`for the normalized`):
print(`r.v. "number of  inversions" of k 1's and n-k 2's up to order s`):
print(`in (1/k), for fixed k and n very large`):
print(`in terms of x`):
print(`AsyPhiW2kFixed(k,2,x);`):

elif nops([args])=1 and op(1,[args])=AsyPhiW2t then
print(`AsyPhiW2t(a,b,t,2,x): the asympt. to order s of the`):
print(`continuous limit `):
print(`for the normalized`):
print(`r.v. "number of  inversions" of a*t 1's and b*t 2's up to order s`):
print(`in terms of. For example, try:`):
print(`AsyPhiW2t(a,b,t,2,x);`):


elif nops([args])=1 and op(1,[args])=AsyRthMomW2t then
print(`AsyRthMomW2t(r,a,b,t,s): the first s terms in the`):
print(`asymptotic expansion, w.r.t. t, of the 2r-th moment `):
print(`(about the mean) of the`):
print(`r.v. Number of inversions for words with a*t 1's and b*t 2's`):
print(`Here a,b,t, and r are symbolic(!), only s is numeric.`):
print(`For example, try:`):
print(`AsyRthMomW2t(r,a,b,t,3); `):


elif nops([args])=1 and op(1,[args])=BetterAppxW2 then
print(`BetterAppxW2(a1,b1,x,s): An approximtion for`):
print(`the number of words in {1^a1,2^b1}  with x inversions.`):
print(`using s terms in the appx. to the continuous limit`):
print(`For example, try: BetterAppxW2(20,20,200,2);`):

elif nops([args])=1 and op(1,[args])=BetterAppxW2kFixed then
print(`BetterAppxW2kFixed(a1,b1,x,s): An approximtion for`):
print(`the number of words in {1^k,2^(n-k)}  with x inversions.`):
print(`using s terms where k is small (fixed) and n is large.`);
print(`For example, try: `):
print(`BetterAppxW2kFixed(50,4,100,2);`):

elif nops([args])=1 and op(1,[args])=Bdok then
print(`Bdok(a1,b1,x,s): compares the allx with one of them small`):
print(`and both big. For example, try:`):
print(`Bdok(200,4,400,2) ;`):

elif nops([args])=1 and op(1,[args])=CheckAppxW2 then
print(`CheckAppxW2(a1,b1):Checks  AppxW2(a1,b1)`):
print(`For example, try: CheckAppxW2(20,20);`):

elif nops([args])=1 and op(1,[args])=CheckBetterAppxW2 then
print(`CheckBetterAppxW2(a1,b1,s):Checks  BetterAppxW2(a1,b1,s)`):
print(`For example, try: CheckBetterAppxW2(20,20,2);`):


elif nops([args])=1 and op(1,[args])=CheckBetterAppxW2kFixed then
print(`CheckBetterAppxW2kFixed(k1,s,T,N1,N2):`):
print(`checks BetterAppxW2kFixed(a1,b1,x,s);`):
print(`for b1=k1 and a1 between N1 and N2`):
print(`CheckBetterAppxW2kFixed(4,2,0,100,110);`):

elif nops([args])=1 and op(1,[args])=CheckithBMk1 then
print(`CheckithBMk1(r,k,K): checks ithBMk(r,a,k) for a[i]'s from 1 to K`):
print(`Do: CheckithBMk1(4,2,10);`):

elif nops([args])=1 and op(1,[args])=CompareAppxWk then
print(`CompareAppxWk(L,x): compares the appx. given by`):
print(`AppxWk(L,x) to the true value, followed by the `):
print(`relative error. For example, try:`):
print(`CompareAppxWk([10,10,5,8],E1(2,[10,10,5,8],4)/2-10);`):

elif nops([args])=1 and op(1,[args])=Fabq then
print(`Fabq(a,b,q): qbinomial(a+b,a)/q^(a*b/2)/(binomial(a+b,a)`):


elif nops([args])=1 and op(1,[args])=FixedkAlpha then
print(`FixedkAlpha(i,k,s): the normalized 2i-th moment`):
print(`also divided by (2*i)!/i!/2^i`):
print(`for the number of inversions with k 1's and n-k 2's`):
print(`with a fixed k, as n goes to infinity, up to s terms`):
print(`in 1/k. For example, try:`):
print(`FixedkAlpha(i,k,s); `):

elif nops([args])=1 and op(1,[args])=Gabq then
print(`Gabq(a,b,q): qbinomial(a+b,a)`):


elif nops([args])=1 and op(1,[args])=ithBM2 then
print(`ithBM2(r,a,b): the polynomial expression,`):
print(`that describes  E[binomial(X-av,r)] where X is the`):
print(`r.v. "number of inversions" on the set of words with`):
print(`a 1's and b 2's`):
print(`and av is the average (ab/2). For example, try:`):
print(`ithBM2(2,a,b);`):

elif nops([args])=1 and op(1,[args])=ithBM3 then
print(`ithBM3(r,a,b,c): the polynomial expression, `):
print(`that describes  E[binomial(X-av,r)] where X is the`):
print(`r.v. "number of inversions" on the set of words with`):
print(`a 1's and b 2's and c 3's`):
print(`and av is the average (ab+bc+ac)/2. For example, try:`):
print(`ithBM3(2,a,b,c);`):

elif nops([args])=1 and op(1,[args])=ithBMk then
print(`ithBMk(r,a,k): Given numeric pos. integers r and k`):
print(`and a symbol r, outputs`):
print(`the polynomial expression, in a[1], ... , a[k]`):
print(`that describes  E[binomial(X-av,r)] where X is the`):
print(`r.v. "number of inversions" on the set of words with`):
print(`a[1] 1's and a[2] 2's ... a[k] k's`):
print(`and av is the average (e_2(a[1], ..., a[k]) /2). For example, try:`):
print(`ithBMk(2,a,3);`):

elif nops([args])=1 and op(1,[args])=ithMomW2 then
print(`ithMomW2(i,a,b): a polynomial expression in a,b, for the i-th moment`):
print(`about the mean for the number of inversions of words `):
print(`with a 1's and b 2's `):
print(`where a and b are symbolic. `):
print(`For example, try:`):
print(`ithMomW2(4,a,b);`):

elif nops([args])=1 and op(1,[args])=ithMomWk then
print(`ithMomWk(i,a,k): a polynomial expression in a[1],..., a[k]`):
print(` for the i-th moment`):
print(`about the mean for the number of inversions of words `):
print(`with a[1] 1's ... a[k]'s k's `):
print(`where a[1], ... a[k] are symbolic. `):
print(`For example, try:`):
print(`ithMomWk(4,a,2);`):

elif nops([args])=1 and op(1,[args])=ithMomW2t then
print(`ithMomW2t(i,a,b,t): ithMomW2(i,a*t,b*t)`):
print(`ithMomW2t(4,a,b,t);`):

elif nops([args])=1 and op(1,[args])=ithMomWktE then
print(`ithMomWktE(i,e,t): the ith-moment about the mean of the number`):
print(`of inversions of a[1]*t 1's ... a[k]*t k's  in terms of the`):
print(`elementary symmetric functions in a[1], ..., a[k].`):
print(`For example, try:`):
print(`ithMomWktE(2,e,t);`):

elif nops([args])=1 and op(1,[args])=ithMomWktM then
print(`ithMomWktM(i,m,t): the ith-moment about the mean of the number`):
print(`of inversions of a[1]*t 1's ... a[k]*t k's  in terms of the`):
print(`monomial symmetric functions m in a[1], ..., a[k].`):
print(`For example, try:`):
print(`ithMomWktM(2,m,t);`):

elif nops([args])=1 and op(1,[args])=ithMomWkt then
print(`ithMomWkt(i,a,k,t): ithMomWk(i,a,k)`):
print(`with a[i] replaced by a[i]*t, for i=1..k`):
print(`For example, try:`):
print(`ithMomWkt(4,a,3,t);`):


elif nops([args])=1 and op(1,[args])=NithMomW2t then
print(`NithMomW2t(i,a,b,t): normalized ithMomW2(i,a*t,b*t)`):
print(`NithMomW2t(4,a,b,t);`):

elif nops([args])=1 and op(1,[args])=NithMomWkt then
print(`NithMomWkt(i,e,t): normalized ith-Moment in terms of`):
print(`elementary symmetric functions`):
print(`NithMomWkt(2,e,t);`):

elif nops([args])=1 and op(1,[args])=PeleW2 then
print(`PeleW2(z,a,b,R):The first R terms of the formal power series`):
print(`((1+z)^-(b/2)-(1+z)^(a+b/2))*a/(1-(1+z)^a/(a+b);`):
print(`For example, try:a`):
print(`PeleW2(z,a,b,5);`):
print(``):

elif nops([args])=1 and op(1,[args])=PeleWk then
print(`PeleWk(z,a,R,k):The first R terms of the formal power series`):
print(`((1+z)^-((a[2]+ ...+a[k]/2)`):
print(`-(1+z)^(a[1]+(a[2]+...+a[k]/2))*a/(1-(1+z)^a[1]/(a[1]+...+a[k]);`):
print(`For example, try:a`):
print(`PeleWk(z,a,5,3);`):
print(``):

elif nops([args])=1 and op(1,[args])=PolToSym then
print(`PolToSym(P,a,k,m): Given a polynomial, P, symmetric`):
print(`in a[1], ..., a[k] and a letter m, condenses it`):
print(`using the monomial symmetric functions m[lamba].`):
print(`For example, try:`):
print(`PolToSym((a[1]+a[2]+a[3]+a[4])^2,a,4,m);`):

else
print(`There is no ezra for`,args):
fi:
 

end:

#PeleW2(z,a,b,R):The first R terms of the formal power series
#((1+z)^-(b/2)-(1+z)^(a+b/2))*a/(1-(1+z)^a/(a+b);
#For example, try:
#PeleW2(z,a,b,5);
PeleW2:=proc(z,a,b,R) local gu,i:
gu:=(a/(a+b))*((1+z)^(-b/2)-(1+z)^(a+b/2))/(1-(1+z)^(a));
gu:=taylor(gu,z=0,R+3):
add(normal(coeff(gu,z,i))*z^i,i=0..R):
end:



#ithBM2(r,a,b): the polynomial expression,
#that describes  E[binomial(X-av,r)] where X is the
#r.v. "number of inversions" on the set of words with
#a 1's and b 2's
#and av is the average (ab/2). For example, try:
#ithBM2(2,a,b);
ithBM2:=proc(r,a,b)
local gu,B,s,z,m:
option remember:
if r=0 then
 RETURN(1):
fi:

gu:=PeleW2(z,a,b,r):
B:=add(coeff(gu,z,s)*subs(a=a-1,ithBM2(r-s,a,b)),s=1..r);
B:=sum(B,a=1..m):
factor(subs(m=a,B)):
end:

Fabq:=proc(a,b,q)
local  i:
expand(normal(mul((1-q^(i+a))/(1-q^i),i=1..b))/q^(a*b/2)/(a+b)!*a!*b!):
end:

Gabq:=proc(a,b,q)
local  i:
expand(normal(mul((1-q^(i+a))/(1-q^i),i=1..b))):
end:



#CheckithBM2(r,K): checks ithBM2(r,a,b) for i and j from1 to K
#Do: CheckithBM2(4,10);
CheckithBM2:=proc(r,K) local gu,i,j,a,b,q:
gu:=ithBM2(r,a,b):

{
seq(
seq(
subs({a=i,b=j}, gu)
 -
subs(q=1,diff(Fabq(i,j,q),q$r)/r!),
i=1..K),j=1..K)}:

end:

ithMomW2:=proc(i,a,b) local gu,r:
gu:=PowerToBinom(i):
factor(add(gu[r+1]*ithBM2(r,a,b),r=0..i)):
end:


ithMomW2t:=proc(i,a,b,t) local gu,i1:
option remember:
gu:=ithMomW2(i,a,b):
gu:=expand(subs({a=t*a,b=t*b},gu)):
add(factor(coeff(gu,t,3*i-i1))*t^(3*i-i1),i1=0..3*i):
end:

#NithMomW2t(i,a,b,t): normalized ith-moment
NithMomW2t:=proc(i,a,b,t) local gu,i1:
option remember:

gu:=ithMomW2t(2*i,a,b,t):
gu:=gu/(a*b*(a+b)/12)^i/((2*i)!/i!/2^i):
add(factor(coeff(gu,t,3*i-i1))*t^(-i1),i1=0..3*i):
end:


#AppxW2(a1,b1,x): An approximtion for
#the number of words in {1^a1,2^b1}  with x inversions.
#For example, try: AppxW2(20,20,200);
AppxW2:=proc(a1,b1,x) local y,sig,a,b:
sig:=evalf(sqrt(subs({a=a1,b=b1},ithMomW2(2,a,b)))):
y:=evalf(a1*b1/2-x)/sig:
if abs(y)>4 then
 RETURN(FAIL):
fi:

evalf((a1+b1)!/a1!/b1!*exp(-y^2/2)/sqrt(2*Pi)/sig):
end:

CheckAppxW2:=proc(a1,b1) local x,av,gu,q,sig,a,b,i:
av:=a1*b1/2:
sig:=evalf(sqrt(subs({a=a1,b=b1},ithMomW2(2,a,b)))):
gu:=expand(normal(mul((1-q^(i+a1))/(1-q^i),i=1..b1))):

{seq(evalf(AppxW2(a1,b1,x)/coeff(gu,q,x)),x=trunc(av-2*sig)..trunc(av+2*sig))}:
end:




#AsyNithMomW2t(r,a,b,t,s): the first s terms in the
#asymptotic expansion, w.r.t. t, of the NORMALIZED 2r-th moment 
#(about the mean) of the
#r.v. Number of inversions for words with a*t 1's and b*t 2's
#Here a,b,t, and r are symbolic(!), only s is numeric.
#For example, try:
#AsyNithMomW2t(a,b,t,r,s)
AsyNithMomW2t:=proc(r,a,b,t,s) local gu,i1,s1,lu,r1,lu1,mu:
gu:=1:
for s1 from 1 to s do
lu:=
expand(
[seq( ((a+b)*a*b)^s1*coeff(NithMomW2t(r1,a,b,t),t,-s1),
r1=1..2*s1+2)]):
mu:=0:
for i1 from 0 to 2*s1 do

lu1:=[seq(coeff(coeff(lu[r1],a,i1),b,2*s1-i1),r1=1..nops(lu))]:

lu1:=factor(GuessPol(lu1,r)):

if lu1=FAIL then
  RETURN(FAIL):
fi:
mu:=mu+lu1*a^i1*b^(2*s1-i1):
od:
gu:=gu+mu*t^(-s1)/((a+b)*a*b)^s1:

od:

for r1 from 1 to 2*s+4  do
if
{seq(normal(coeff(NithMomW2t(r1,a,b,t),t,-s1)
-coeff(subs(r=r1,gu),t,-s1)),s1=1..s)}<>{0} then

print(`Something is wrong`):
RETURN(FAIL):
fi:
od:

gu:
end:


#AsyRthMomW2t(R,a,b,t,s): the first s terms in the
#asympt. expression in t for the 2r-th moment about
#the mean for the r.v. NumberOfIntversions for words
#in t*a, t*b
AsyRthMomW2t:=proc(r,a,b,t,s)
t^(3*r)*AsyNithMomW2t(r,a,b,t,s)*(2*r)!/r!/2^r*(a*b*(a+b)/12)^r:
end:

#AsyAlphaW2tS(r,a,b,t,s): The first R terms in the
#asymptotic expansion for alpha_{2r} (normlaized 2r-th moment)
#for the r.v. NumberOfInversions on words with a*t 1's and b*t 2's
#For example, try:
#AsyAlphaW2tS(R,a,b,t,2);
AsyAlphaW2tS:=proc(r,a,b,t,s) local gu,gu1,x,i1:
#(2*r)!/r!/2^r:
gu:=normal(subs(t=1/x,AsyNithMomW2t(r,a,b,t,s))):
gu1:=normal(subs(t=1/x,NithMomW2t(1,a,b,t))):
gu:=taylor(gu/gu1^r,x=0,s+1):
gu:=add(factor(coeff(gu,x,i1))/t^i1,i1=0..s):
(2*r)!/r!/2^r*collect(gu,t):
end:



#AlphaW2(r1,a,b,t): alpha_{2r1} of the r.v. Number of Inversions
#with a*t 1's and b*t 2's for numeric r1. For example, try:
#AlphaW2t(3,a,b,t);
AlphaW2t:=proc(r1,a,b,t):
normal(ithMomW2t(2*r1,a,b,t)/ithMomW2t(2,a,b,t)^r1):
end:


#AsyAlphaW2t(r,a,b,t,s): The first s terms in the
#asymptotic expansion for alpha_{2r} (normlaized 2r-th moment)
#for the r.v. NumberOfInversions on words with a*t 1's and b*t 2's
#and numeric r
#For example, try:
#AsyAlphaW2t(2,a,b,t,2);
AsyAlphaW2t:=proc(r1,a,b,t,s) local gu,x,i:
gu:=AlphaW2t(r1,a,b,t):
gu:=normal(subs(t=1/x,gu)):
gu:=taylor(gu,x=0,s+1):
add(normal(coeff(gu,x,i))/t^i,i=0..s):
end:




#The asymptotic moment generating function of the normalized
#r.v. "number of  inversions" of a*t 1's and b*t 2's up to order s
#in terms of x
#AsyMgfW2t(a,b,t,2,x);
AsyMgfW2t:=proc(a,b,t,s,x) local gu,r,i:
gu:=AsyAlphaW2tS(r,a,b,t,s):
add(sum(coeff(gu,t,-i)*x^(2*r)/(2*r)!,r=0..infinity)/t^i,i=0..s):
end:


#AsyPhiW2t(a,b,t,2,x): the asympt. to order s of the
#continuous limit 
#for the normalized
#r.v. "number of  inversions" of a*t 1's and b*t 2's up to order s
#in terms of x
#AsyPhiW2t(a,b,t,2,x);
AsyPhiW2t:=proc(a,b,t,s,x) local gu,X,i1:
gu:=AsyMgfW2t(a,b,t,s,x):
gu:=subs(x=-I*x,gu):
gu:=
add(int(coeff(gu,t,-i1)*exp(I*x*X),x=-infinity..infinity)/t^i1,
i1=0..s)/2/Pi:
subs(X=x,gu):
end:





#BetterAppxW2(a1,b1,x,s): An approximtion for
#the number of words in {1^a1,2^b1}  with x inversions.
#using s terms
#For example, try: 
#BetterAppxW2(20,20,200,2);
BetterAppxW2:=proc(a1,b1,x,s) local y,sig,a,b,X,gu,t:
gu:=AsyPhiW2t(a,b,t,s,X):
gu:=subs({a=a1/(a1+b1),b=b1/(a1+b1),t=a1+b1},gu):
sig:=evalf(sqrt(subs({a=a1,b=b1},ithMomW2(2,a,b)))):
y:=evalf(a1*b1/2-x)/sig:
if abs(y)>4 then
 RETURN(FAIL):
fi:

evalf((a1+b1)!/a1!/b1!*subs(X=y,gu)/sig):

end:

CheckBetteArppxW2:=proc(a1,b1,s) local x,av,gu,q,sig,a,b,i:
av:=a1*b1/2:
sig:=evalf(sqrt(subs({a=a1,b=b1},ithMomW2(2,a,b)))):
gu:=expand(normal(mul((1-q^(i+a1))/(1-q^i),i=1..b1))):

{seq(evalf(
BetterAppxW2(a1,b1,x,s)/coeff(gu,q,x)),x=trunc(av-2*sig)..trunc(av+2*sig))}:
end:


#########start W3
#PeleW3(z,a,b,c,R):The first R terms of the formal power series
#((1+z)^-(b/2)-(1+z)^(a+b/2))*a/(1-(1+z)^a/(a+b);
#For example, try:
#PeleW3(z,a,b,c,5);
PeleW3:=proc(z,a,b,c,R) local gu,i:
gu:=(a/(a+b+c))*((1+z)^(-(b+c)/2)-(1+z)^(a+(b+c)/2))/(1-(1+z)^(a));
gu:=taylor(gu,z=0,R+3):
add(normal(coeff(gu,z,i))*z^i,i=0..R):
end:

#ithBM3(r,a,b,c): the polynomial expression,
#that describes  E[binomial(X-av,r)] where X is the
#r.v. "number of inversions" on the set of words with
#a 1's and b 2's and c 3's
#and av is the average ((ab+ac+bc)/2). For example, try:
#ithBM3(2,a,b,c);
ithBM3:=proc(r,a,b,c)
local gu,B,s,z,m:
option remember:
if r=0 then
 RETURN(1):
fi:

gu:=PeleW3(z,a,b,c,r):
B:=add(coeff(gu,z,s)*subs(a=a-1,ithBM3(r-s,a,b,c)),s=1..r);
B:=sum(B,a=1..m)+ithBM2(r,b,c):
factor(subs(m=a,B)):
end:



#PeleWk(z,a,R,k):The first R terms of the formal power series
#((1+z)^-(e_1(a[2], ...a[k])/2)-(1+z)^(a[1]+(a[2]+...+a[k]/2))
#*a[1]/(1-(1+z)^a[1]/(a[1]+a[2]+ ...);
#For example, try:
#PeleWk(z,a,5,k);
PeleWk:=proc(z,a,R,k) local gu,i,i1:
gu:=(a[1]/
(add(a[i1],i1=1..k)))*
((1+z)^(-(add(a[i1],i1=2..k))/2)-
(1+z)^(a[1]+(add(a[i1],i1=2..k))/2))/(1-(1+z)^(a[1]));
gu:=taylor(gu,z=0,R+3):
add(normal(coeff(gu,z,i))*z^i,i=0..R):
end:

#ithBMk(r,a,k): Given numeric pos. integers r and k
#and a symbol r, outputs
#the polynomial expression, in a[1], ... , a[k]
#that describes  E[binomial(X-av,r)] where X is the
#r.v. "number of inversions" on the set of words with
#a[1] 1's and a[2] 2's ... a[k] k's
#and av is the average (e_2(a[1], ..., a[k]) /2). For example, try:
#ithBMk(2,a,3);
ithBMk:=proc(r,a,k)
local gu,B,s,z,m,i,lu:
option remember:
if r=0 then
 RETURN(1):
fi:

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


gu:=PeleWk(z,a,r,k):
B:=add(coeff(gu,z,s)*subs(a[1]=a[1]-1,ithBMk(r-s,a,k)),s=1..r);
lu:=ithBMk(r,a,k-1):
lu:=subs({seq(a[i]=a[i+1],i=1..k-1)},lu):
B:=sum(B,a[1]=1..m)+lu:
factor(subs(m=a[1],B)):
end:




Fkq:=proc(a,q)
local  i,j,mu,lu,i1,i2:
mu:=add(add(a[i1]*a[i2],i2=i1+1..nops(a)),i1=1..nops(a)):
mu:=q^(mu/2):
mu:=mu*convert(a,`+`)!/mul(a[i]!,i=1..nops(a)):
lu:=expand(normal(mul(1-q^i,i=1..convert(a,`+`))/
mul(mul(1-q^i,i=1..a[j]),j=1..nops(a)))):
expand(lu/mu):
end:

#CheckithBMk1(r,k,K): checks ithBMk(r,a,k) for a[i]'s from 1 to K
#Do: CheckithBMk1(4,2,10);
CheckithBMk1:=proc(r,k,K) local gu,a,q,ra,lu,i1:
gu:=ithBMk(r,a,k):
ra:=rand(1..K):
lu:=[seq(ra(),i1=1..k)]:

evalb(
subs({seq(a[i1]=lu[i1],i1=1..k)}, gu)
 =
subs(q=1,diff(Fkq(lu,q),q$r)/r!)):
end:


ithMomWk:=proc(i,a,k) local gu,r:
gu:=PowerToBinom(i):
expand(add(gu[r+1]*ithBMk(r,a,k),r=0..i)):
end:


ithMomWkt:=proc(i,a,k,t) local gu,i1:
option remember:
gu:=ithMomWk(i,a,k):
gu:=expand(subs({seq(a[i1]=a[i1]*t,i1=1..k)},gu)):
add(factor(coeff(gu,t,3*i-i1))*t^(3*i-i1),i1=0..3*i):
end:

Can1:=proc(a) local i:
for i from 1 to nops(a) while a[i]=0 do od:
[seq(a[nops(a)-j],j=0..nops(a)-i)]:
end:

#PolToSym(P,a,k,m): Given a polynomial, P, symmetric
#in a[1], ..., a[k] and a letter m, condenses it
#using the monomial symmetric functions m[lamba]
#For example, try:
#PolToSym((a[1]+a[2]+a[3]+a[4])^2,a,4,m);
PolToSym:=proc(P,a,k,m) local P1,M,P11,lu,C,i1,gu,i,lu1:
P1:=expand(P):

if not type(P1,`+`) then
 RETURN(FAIL):
fi:

M:={}:
gu:=0:

for i from 1 to nops(P1) do
P11:=op(i,P1):
lu:=[seq(degree(P11,a[i1]),i1=1..k)]:
lu1:=sort(lu):
 if not member(lu1,M) then
    M:=M union {lu1}:
 C:=normal(P11/mul(a[i1]^lu[i1],i1=1..k)):
  if not type(C,numeric) then
     RETURN(FAIL):
  fi:
  gu:=gu+C*m[op(Can1(lu1))]:
 fi:
od:
gu:
end:

    




#ithMomWktM(i,m,t): the ith-moment about the mean of the number
#of inversions of a[1]*t 1's ... a[k]*t k's  in terms of the
#monomial symmetric functions m in a[1], ..., a[k].
#For example, try:
#ithMomWktM(2,m,t);
ithMomWktM:=proc(i,m,t) local k,a,d,gu,lu,i1,vu:
if i mod 2 =1 then
 RETURN(0):
fi:

k:=3*i/2+1:

gu:=expand(ithMomWkt(i,a,k,t)):


d:=degree(gu,t):

vu:=0:
for i1 from 0 to d do
lu:=coeff(gu,t,i1):
if lu<>0 then
lu:=PolToSym(lu,a,k,m):
vu:=vu+lu*t^i1:
fi:
od:

vu:
end:



###begin symmetric functions stuff

#ez:=proc(): print(`mtoe(lam,e), E1(i,x,n), Mekadem(F,x,lam)`): end:


#E1(i,x,n): the ith-elementary symmetric function in n variables
E1:=proc(i,x,n) local t,j:
coeff(expand(mul(1+t*x[j],j=1..n)),t,i):
end:


#Mekadem(F,x,lam): the coeff. of x[1]^lam[1]*x[2]^lam[2]...
Mekadem:=proc(F,x,lam) local gu,i:

gu:=F:

for i from 1 to nops(lam) do
gu:=coeff(gu,x[i],lam[i]):
od:
gu:
end:






#mtoe(lam,e): expresses m[lam] as a polynomial in the elementary
#symmetric functions e[i]'s
mtoe:=proc(lam,e) local n,mu,gu,a,eq,x,i,j,lam1,gu1,var:
n:=convert(lam,`+`):
mu:=partition(n):
gu:=0:
gu1:=0:
var:={}:
for i from 1 to nops(mu) do
var:=var union {a[i]}:
lam1:=mu[i]:
gu:=gu+a[i]*mul(e[lam1[j]],j=1..nops(lam1)):
gu1:=expand(gu1+a[i]*mul(E1(lam1[j],x,n),j=1..nops(lam1))):
od:

eq:={}:

for i from 1 to nops(mu) do
lam1:=mu[i]:
 if lam1=lam then
    eq:=eq union {Mekadem(gu1,x,lam1)-1}:
 else
   eq:=eq union {Mekadem(gu1,x,lam1)}:
  fi:
od:

var:=solve(eq,var):

if var=FAIL then
 RETURN(FAIL):
fi:

subs(var,gu):

end:

###end symmetric functions stuff


Rev:=proc(L) local i: [seq(L[nops(L)+1-i],i=1..nops(L))]:end:

#ithMomWktE(i,e,t): the ith-moment about the mean of the number
#of inversions of a[1]*t 1's ... a[k]*t k's  in terms of the
#elementary symmetric functions m in a[1], ..., a[k].
#For example, try:
#ithMomWktE(2,e,t);
ithMomWktE:=proc(i,e,t) local gu,m,lu,mu,j,lu1,lam,lu11,i1:
if i mod 2 =1 then
 RETURN(0):
fi:
lu:=ithMomWktM(i,m,t):
gu:=0:

for j from 1 to degree(lu,t) do

lu1:=coeff(lu,t,j):

mu:=partition(j):

for lam in mu do
 lu11:=coeff(lu1,m[op(Rev(lam))],1):
 gu:=gu+lu11*mtoe(lam,e)*t^j:
od:

od:

add(factor(coeff(gu,t,degree(gu,t)-i1))*t^(degree(gu,t)-i1),i1=0..degree(gu,t))
:
end:





#NithMomWkt(i,e,t): normalized ith-moment for the multivariable case
#in terms of the elementary symmetric functions
NithMomWkt:=proc(i,e,t) local gu,i1,gu1:
option remember:

gu1:=coeff(ithMomWktE(2,e,t),t,3):
gu:=ithMomWktE(2*i,e,t):

gu:=gu/(gu1)^i/((2*i)!/i!/2^i):
add(factor(coeff(gu,t,3*i-i1))*t^(-i1),i1=0..3*i):
end:






#FixedkAlpha(i,k,s): the normalized 2i-th moment
#for the number of inversions with k 1's and n-k 2's
#with a fixed k, as n goes to infinity, up to s terms
#in 1/k. For example, try:
#FixedkAlpha(i,k,s) 
FixedkAlpha:=proc(i,k,s) local gu,a,b,gu1,i1,n:
gu:=coeff(subs({a=k,b=n-k},ithMomW2(2*i,a,b)),n,2*i):
gu1:=coeff(subs({a=k,b=n-k},ithMomW2(2,a,b)),n,2):
gu:=gu/gu1^i:
gu:=gu/((2*i)!/i!/2^i):
add(coeff(gu,k,-i1)/k^i1,i1=0..s):
end:





#AsyFixedkAlpha(r,k,s): the asymptotic normalized 2i-th moment
#for the number of inversions with k 1's and n-k 2's
#with a fixed k, as n goes to infinity, up to s terms
#in 1/k. For example, try:
#AsyFixedkAlpha(r,k,2);
AsyFixedkAlpha:=proc(r,k,s) local gu,i,j:

gu:=[seq(FixedkAlpha(i,k,s),i=1..2*s+2)]:
add(factor(GuessPol([seq(coeff(gu[i],k,-j),i=1..2*s+2)],r))/k^j,j=0..s):
end:



#AsyMgfW2kFixed(k,x,s):
#The asymptotic moment generating function of the normalized
#r.v. "number of  inversions" of k 1's and n-k 2's up to order s
#with k fixed and n large
#in terms of x. For example, try:
#AsyMgfW2kFixed(k,x,2);
AsyMgfW2kFixed:=proc(k,x,s) local gu,r,i:
gu:=AsyFixedkAlpha(r,k,s):
add(sum(coeff(gu,k,-i)*x^(2*r)/r!/2^r
,r=0..infinity)/k^i,i=0..s):
end:



#AsyPhiW2kFixed(k,2,x): the asympt. to order s of the
#continuous limit 
#for the normalized
#r.v. "number of  inversions" of k 1's and n-k 2's up to order s
#in (1/k), for fixed k and n very large
#in terms of x
#AsyPhiW2kFixed(k,2,x);
AsyPhiW2kFixed:=proc(k,s,x) local gu,X,i1:
gu:=AsyMgfW2kFixed(k,x,s):
gu:=subs(x=-I*x,gu):
gu:=
add(int(coeff(gu,k,-i1)*exp(I*x*X),x=-infinity..infinity)/k^i1,
i1=0..s)/2/Pi:
subs(X=x,gu):
end:



#BetterAppxW2kFixed(a1,b1,x,s): An approximtion for
#the number of words in {1^k,2^(n-k)}  with x inversions.
#using s terms where k is small (fixed) and n is large
#For example, try: 
#BetterAppxW2kFixed(50,4,100,2);
BetterAppxW2kFixed:=proc(a1,b1,x,s) local y,sig,X,gu,k,a,b,k1:
k1:=min(a1,b1):


gu:=subs(k=k1,AsyPhiW2kFixed(k,s,X)):

sig:=evalf(sqrt(subs({a=a1,b=b1},ithMomW2(2,a,b)))):
y:=evalf(a1*b1/2-x)/sig:

if abs(y)>4 then
 RETURN(FAIL):
fi:


evalf((a1+b1)!/a1!/b1!*subs(X=y,gu)/sig):

end:




#CheckBetterAppxW2kFixed(k1,s,T,N1,N2):
#checks BetterAppxW2kFixed(a1,b1,x,s);
#for b1=k1 and a1 between N1 and N2
#CheckBetterAppxW2kFixed(4,2,0,100,110);
CheckBetterAppxW2kFixed:=proc(k,s,T,N1,N2) local n1,q,lu:


lu:=seq(
coeff(Gabq(k,2*n1,q),q,k*n1-T)/BetterAppxW2kFixed(k,2*n1,k*n1-T,s),
n1=N1..N2):

min(lu),max(lu):

end:



#Bdok(a1,b1,x,s): compares the allx with one of them small
#and both big. For example, try:
#Bdok(200,4,400,2) ;
Bdok:=proc(a1,b1,x,s) local gu,q:
gu:=coeff(Gabq(a1,b1,q),q,x):
[(BetterAppxW2kFixed(a1,b1,x,s)-gu)/gu,
(BetterAppxW2(a1,b1,x,s)-gu)/gu,gu]:
end:








#AlphaWktE(r1,e,t): alpha_{2r1} of the r.v. Number of Inversions
#with a[1]*t 1's, a[2] 2's ..., a[n]*t n's for numeric r1. For example, try:
#AlphaWktE(3,e,t);
AlphaWktE:=proc(r1,e,t):
normal(ithMomWktE(2*r1,e,t)/ithMomWktE(2,e,t)^r1):
end:


#AsyAlphaWktE(r,e,t,s): The first s terms in the
#asymptotic expansion for alpha_{2r} (normlaized 2r-th moment)
#for the r.v. NumberOfInversions on words with a[1]*t 1's 
#a[2]*t 2's ... a[n] n's in terms of the elementary symmtric
#functions in a[1], ..., a[n], denoted by e[i]'s
#r and s are numeric, t and e are symbolic
#For example, try:
#AsyAlphaWktE(2,e,t,2);
AsyAlphaWktE:=proc(r1,e,t,s) local gu,x,i:
gu:=AlphaWktE(r1,e,t):
gu:=normal(subs(t=1/x,gu)):
gu:=taylor(gu,x=0,s+1):
add(normal(coeff(gu,x,i))/t^i,i=0..s):
end:


#AppxWk(L,x): An approximtion for
#the number of words in {1^L[1],2^L[2], ..., n^L[n]}  with x inversions.
#For example, try: AppxWk([20,20,20],600);
AppxWk:=proc(L,x) local y,sig,gu,n,t,L1,t1,e,i:
n:=nops(L):
gu:=ithMomWktE(2,e,t);
t1:=convert(L,`+`):
L1:=[seq(L[i]/t1,i=1..nops(L))]:
gu:=subs({t=t1,e[1]=E1(1,L1,n),e[2]=E1(2,L1,n),e[3]=E1(3,L1,n)},gu):
sig:=evalf(sqrt(gu)):
y:=evalf(E1(2,L,n)/2-x)/sig:
if abs(y)>4 then
 RETURN(FAIL):
fi:

evalf(convert(L,`+`)!/mul(L[i]!,i=1..n)*exp(-y^2/2)/sqrt(2*Pi)/sig):
end:


#CompareAppxWk(L,x): compares the appx. given by
#AppxWk(L,x) to the true value, followed by the 
#relative error. For example, try:
#CompareAppxWk([10,10,5,8],E1(2,[10,10,5,8],4)/2-10);
CompareAppxWk:=proc(L,x) local gu,n,q,gu1,i,j:
n:=nops(L):
gu:=normal(
mul(1-q^i,i=1..convert(L,`+`))/mul(mul(1-q^j,j=1..L[i]),i=1..n)):
gu:=coeff(gu,q,x):
gu1:=AppxWk(L,x):
gu,gu1,abs(gu-gu1)/gu:
end: