# restart; read `c:/Users/Asus/Google Drive
# /Aek/MomentReview/InvMaj.txt`;
# read `c:/Users/Asus/Google Drive/Guess.txt`;

with(combinat,permute):
with(combinat,stirling2):

#print(`Welcome to InvMaj, a maple package `):
#print(`written by Thotsaporn (Aek) Thanatipanonda`):
#print(`First Version: September 24, 2013 `):
#print(`Second Version: March 3, 2020 `):
#print(``):
#print(`It is a simpler ananlog of`): 
#print(`Dr.Z with Baxter,`):
#print(`The Number of Inversions and `);
#print(`the Major Index of Permutations are `):
#print(`Asymptotically Joint-Independently Normal`):

Help := proc()

print(`Objects: Inv and Maj`);
print(`Section 1: Empirical`);
print(`The functions are`):
print(`Inv(S), Maj(S), `);
print(`Moment(r,s,n), MomentMean(r,s,n),`);
print(`ForMo(r,s,n), ForMoM(r,s,n),`);
print(`Corre(n)`);
print(``);
print(`Section 2: Asymp. behavior of inv`);
print(`InvFn(n,q), InvGn(n,q), InvB(r,n)`);
print(`Pn(n,q), InvPs(r,n), `);
print(`FormuInvPs(r,n), InvLeadB(r,n,t),`);
print(`RecB(r,n), MoInv(r,n),`);
print(`InvLeadM(r,n,t), ZMoInv(m,r,n)`);
print(``);
print(`Section 3: Asymp. behavior of maj`);
print(`MajF(n,i,q), MajH(n,q), MajG(n,i,q),`);
print(`MajBG(R,n,i), GuessMajB(r,n,i),`);
print(`TestMaj4(n,i,q), CofEx4(n),`);
print(`MajVer6(r)`);
print(``);

print(`For joint distribution of inv and maj,`);
print(`look at program Baxter.`);
print(``);

end:




########################
# Objects: Inv and Maj

# Section 1: Empirical
# Input: permutation of word.
# Ouput: number of inversion of sigma.
# Try: Inv([3,4,2,1]);
Inv := proc(S) option remember;
local i,n,count;

count := 0;
n := nops(S);

if n =1 then
	return(0);
fi:

for i from 1 to n-1 do
	if S[i] > S[n] then
		count := count+1;
	fi:
od:

count+Inv([op(1..n-1,S)]);
end:


# Input: permutation sigma.
# Ouput: number of major index of sigma.
# Try: Maj([3,4,2,1]);
Maj := proc(S) option remember;
local i,n,maj;

n := nops(S);
maj := 0;

for i from 1 to n-1 do
	if S[i]>S[i+1] then
		maj := maj+i;
	fi:
od:

maj;
end:


####################################
# Average on perm. of n
# of inv^r*maj^s by brute force.
# Try: seq(Moment(1,0,n),n=1..5);
Moment := proc(r,s,n) option remember;
local p;
add(Inv(p)^r*Maj(p)^s
, p in permute(n))/n!;
end:


# Average on perm. of n
# of (inv-E(inv))^r*(maj-E(maj))^s
# by brute force.
# Try: seq(MomentMean(1,1,n),n=1..5);
MomentMean := proc(r,s,n) option remember;
local p,minv,mmaj;

minv := Moment(1,0,n);
mmaj := Moment(0,1,n);
add((Inv(p)-minv)^r*(Maj(p)-mmaj)^s,
p in permute(n))/n!;

end:


# Compute formula in n, the average 
# of inv^r*maj^s 
# over all perm(n) by brute force.
# Try: ForMo(1,1,n);
ForMo := proc(r,s,n) option remember;
local A,i;

A := [seq(Moment(r,s,i),i=1..2*(r+s)+1)];
FitPoly(A,1,n);
end:


# Compute formula in n, the average 
# of (inv-E(inv))^r*(maj-E(maj))^s 
# over all perm(n) by brute force.
# Try: ForMoM(1,1,n);
# seq(ForMoM(r,0,n),r=1..7);
ForMoM := proc(r,s,n) option remember;
local A,i;

A := [seq(MomentMean(r,s,i),i=1..3/2*(r+s)+2)];
FitPoly(A,1,n);
end:


##########################
# Some interesting Outputs

# Correlation of inv and maj.
# Try: Corre(n);
Corre := proc(n) option remember;

simplify(ForMoM(1,1,n)/ForMoM(2,0,n)^(1/2)
/ForMoM(0,2,n)^(1/2));
end:



#########################
# Section 2: 
# Probability dist. of # 
# of inv. is asymptotically 
# normal.

# Probability Generating Function.
# Try: InvFn(2,q);
InvFn := proc(n,q) 
local i;
factor(1/n!*product((1-q^i)/(1-q),i=1..n));
end:

#Try: InvGn(7,q);
InvGn := proc(n,q) 
InvFn(n,q)/q^(n*(n-1)/4);
end:

# Input: numeric r and n
# Try: InvB(10,5);
InvB := proc(R,n) option remember;
local T,z,i;

T := simplify(taylor(InvGn(n,1+z),z=0,R+2));  
T := convert(T,polynom); 
[1,seq(factor(coeff(T,z^i)),i=1..R)];
end:

# Try: taylor(Pn(n,1+z),z=0,10);
Pn := proc(n,q)
1/n*(1-q^n)/(1-q)*q^((1-n)/2);
end:


# Output: List of coefficients of P.
#Try: InvPs(4,n);
InvPs := proc(r,n) option remember;
local z,T,P,i:

P := Pn(n,1+z);
T := taylor(P,z=0,r+2); 
T := simplify(convert(T,polynom));
[1,seq(factor(coeff(T,z^i)),i=1..r)];

end:


# Output: List of coefficients of P. by formula
#Try: [seq(FormuInvPs(i,n),i=0..7)];
FormuInvPs := proc(r,n) option remember;
local A,s:

A:=1/n*add((-1)^s*MyBi((n-3)/2+s,s)
*MyBi(n,r+1-s),s=0..r);

factor(A);
end:

# Input: numeric r and symbolic n
# Output: the leading term of B_r(n)
# Try: InvLeadB(5,n,0);
InvLeadB := proc(r,n,t) option remember;
local i,A,deg,S;

A := InvB(r,n);
S := [];

for i from 1 to nops(A) do
	deg := degree(A[i])-t;
	S := [op(S),coeff(A[i],n,deg)*n^deg];
od:

S;
end:




# B_r(n) from recursion identity
# Try: RecB(5,n); InvB(5,n);
# t:=time():mm:=[seq(RecB(r,n),r=0..14)]:time()-t;
RecB := proc(r,n) option remember:
local s,P1,DE;	

if r = 0 then
	return(1);
elif r = 1 then
	return(0);
fi:

P1 := InvPs(r,n);
DE := add(subs(n=n-1,RecB(r-s,n))*P1[s+1],s=2..r);

factor(Del(DE,n));
end:



# Moment about the mean of number of inv. 
# by generating function.
# Try: seq(MoInv(r,n),r=1..5);
MoInv := proc(r,n) option remember;
local i,SB;

if r = 0 then
	return(1);
fi:

SB := InvB(r,n);
factor(add(stirling2(r,i)*SB[i+1]*i!,i=0..r));
end:


# Input: numeric r and symbolic n
# Output: the leading term of B_r(n)
# Try: InvLeadM(5,n,0);
InvLeadM := proc(r,n,t) option remember;
local i,A,deg,S;

A := [seq(MoInv(i,n),i=0..r)];
S := [];

for i from 1 to nops(A) do
	deg := degree(A[i])-t;
	S := [op(S),coeff(A[i],n,deg)*n^deg];
od:

S;
end:


# The moment of inversion in r and n
# which accurate up to an (m+1)_th order.
# Try: ZMoInv(3,r,n);
ZMoInv := proc(m,r,n)
local i,s,T,C,sol;

C := (i,n) -> MoInv(2,n)^i*(2*i)!/(2^i*(i!));
sol := 0;

for s from 0 to m do	
	T := [seq(
	GetCoef(MoInv(2*i,n),C(i,n),-s,0,n)
	,i=0..2*s+2)]; 
	sol := sol+GuessPol(T,0,r)/n^s;
od:

C(r,n)*sol;
end:


#########################
# Section 3: 
# Probability dist. of # 
# of maj. is asymptotically 
# normal.

# Output: F(n,i)(q), the generating
# function of major index where the 
# permutation ended in i
# Try: MajF(6,6,q);
MajF := proc(n,i,q) option remember;
local j,A;

if n =1 then
	A := 1;
elif i = n then
	A :=add(MajF(n-1,j,q),j=1..n-1);
else
	A := MajF(n,i+1,q)
	-(1-q^(n-1))*MajF(n-1,i,q);
fi:

factor(A);
end:


# Output: H(n)(q), the generating
# function of major index 
# Try: [seq(MajH(i,q),i=1..5)];
MajH := proc(n,q) option remember;
sort(factor(MajF(n+1,n+1,q)));
end:


# Output: G(n,i)(q), the prob. generating
# function of major index about the
# mean where the permutation ended in i
# Note: G(n,i)(q) are the same for all i.
# Try: 
# n:=7:   MajG(n,n,q);
# 1/(n-1)*factor(add(q^(n/2-j)*MajG(n-1,j,q),j=1..n-1));
# i:=5:  MajG(n,i,q);
# factor( MajG(n,i+1,q)/q
# +(q^(n-1)-1)/q^(n/2)/(n-1)*MajG(n-1,i,q));
MajG := proc(n,i,q) option remember;
MajF(n,i,q)/(n-1)!/q^(n-i+(n-1)*(n-2)/4);

end:


# Output: List of coefficients of G,
# the Binomial Moment.
# Try: MajBG(4,5,3);
# InvB(14,5); MajBG(14,6,6);
MajBG := proc(R,n,i) option remember;
local j,z,T,P:

P := MajG(n,i,1+z);
T := taylor(P,z=0,R+2); 
T := simplify(convert(T,polynom));
[1,seq(factor(coeff(T,z^j)),j=1..R)];
end:


# Output: Guess formula for B_r(n,i) 
# Note: B_r(n,i) is the same for all i. 
# Try: GuessMajB(5,n,2);  InvB(5,n-1);
GuessMajB := proc(r,n,i) option remember;
local j,A;

A :=[seq(MajBG(r+1,j,i)[r+1],j=i..i+20)];
GuessPol(A,i,n);
end:


# Test identity 4 in my note!
# Try: TestMaj4(5,3,q);
TestMaj4 := proc(n,i,q) option remember;
local L,R;

if i >= n then
	ERROR("BadInput");
fi:

L := MajG(n,i,q);
R := 1/q*MajG(n,i+1,q)+(q^(n-1)-1)
/q^(n/2)/(n-1)*MajG(n-1,i,q);

[L,simplify(L-R)];
end:


# Coefficient of the expandsion in eq4.
# Try: CofEx4(n);
CofEx4 := proc(n) option remember;
local z,A;

A :=taylor( ((1+z)^(n-1)-1)/(1+z)^(n/2)/(n-1) ,z,10);
[seq( factor(coeff(A,z,i)),i=0..10)];
end:


# Verify the correctness of id. 4,6 in the note
# Try:  MajVer6(5); 
MajVer6 := proc(r) option remember;
local i,n,T,L,R;

T := CofEx4(n);
L := GuessMajB(r,n);
R := add( (-1)^i*GuessMajB(r-i,n),i=0..r)+
add( T[i+1]*GuessMajB(r-i,n-1),i=0..r);

[L,factor(L-R)];
end: