#OK to post homework

#Zidong Zhang, 28/2/2021, Assignment 10


#Inv(N):the number of inversion, denoted by inv(¦Ð), is the number of pairs 1 ¡Ü i < j ¡Ü n such ¦Ð[i] > ¦Ð[j][For example inv(631425)=8]


Inv := proc(N) local i , j , n , count ,m , L, k, c, b:
n := convert(N , base , 10):
count := 0:

for i from 1 to nops(n)-1 do 
	for j from i+1 to nops(n)  do
		if n[i]<n[j] then
			count := count +1:
		fi:

	od:
od:
count:
end:

#Des(N):The number of descents, denoted by des(¦Ð), is the number of places i, such that ¦Ð[i] > ¦Ð[i+1] [For example inv(631425)=3]


Des := proc(N) local n , i, count :

n := convert(N, base , 10):
count := 0:

for i from 1 to nops(n)-1 do 
	if n[i] < n[i+1] then
		count := count +1:
	fi:
od:
count:
end:

#Maj(N):The major index, denoted by maj(¦Ð), is the SUM of the places i, such that ¦Ð[i] > ¦Ð[i+1] [For example maj(631425)=1+2+4=7]



Maj := proc(N) local n , su, i:

n:= convert(N , base , 10):
su:=0:

for i from 1 to nops(n)-1 do
	if n[i]<n[i+1] then
		su:= su+(nops(n)-i):
	fi:
od:
su:
end:


#ExactGFinv(n,x):inputs a positive integer n and a variable name n, and outputs the generating function whose coefficient of x^k is the exact number of permutations with inv=k


##### I'm kind confused of the question. My understanding is when we try to have the coeffcient, it is not canonical. i.e. there are a lot of permutation that have inv=k, maj=k, des=k. e.g inv(21345) = inv(12354).
#ExactGFinv := proc(n,x) local  :






#end:











#PnSeq(N,t):That inputs a positive integer N and outputs the first N terms of Pn(t)

PnSeq := proc(N,t)  local T,L,i,p:

T:= taylor(exp(t*(exp(x)-1)), x = 0 , N):
L:=[]:
for i from 0 to N do 
	p := coeff(T,x,i)*i!:
	L:=[op(L), p]:
od:
L:
end:

######after comparing it with SAc(exp(t*(exp(x)-1) ),x,20), and not normalized it, I think they are getting closer and closer.


###############from C10.txt


SAc:=proc(f,x,K) local i,mu,sig2,M,f1:
if subs(x=1,f)=0 then
 RETURN(0):
fi:

f1:=f/subs(x=1,f):
mu:=subs(x=1,diff(f1,x)):
f1:=f1/x^mu:
f1:=x*diff(f1,x):
f1:=x*diff(f1,x):
sig2:=subs(x=1,f1):

M:=[mu,sig2]:

if sig2=0 then RETURN(M): fi:
for i from 3 to K do
  f1:=x*diff(f1,x):
#  M:=[op(M), evalf(subs(x=1,f1)/sig2^(i/2))]:
  M:=[op(M), subs(x=1,f1)/sig2^(i/2)]:
od:
M:
end: