#M22.txt: Maple Code for Dr. Z.'s Lecture 22 of Combinatorics, taught at Rutgers University

Help22:=proc(): print(` TreeSeq(N), ATreeSeq(N,r), FunEqToSeq(PHI,z,N), FunEqToSeq(PHI,z,N) , LIF(PHI,z,N), TreeSeq1(N) , TreeSeqL(N,t) `):
print(`AveAndMoms(f,x,K);`):
 end:


#TreeSeq(N): The first N terms of the sequence eunerating the number of labelled trees on n vertices
#alias: "connected labelled graphs in n vertices and n-1 edges"
#Using the edge-weighted egf for labelled connected graphs
#log(Sum((1+a)^binomial(i,2)*x^i/i!,i=0..infinity);
TreeSeq:=proc(N) local f,x,a,n:

#Since we are only interested in the first N terms, we only need the terms up to x^(N) (the rest is irrelenat to us)
f:=log(add((1+a)^(n*(n-1)/2)*x^n/n!,n=0..N+1)):

#we take the Taylor exapansion up to x^N (the rest is garbage)
f:=taylor(f,x=0,N+1):

#We take the coefficient of x^n, multiply by n! (since it is an EXPONENTIAL generating function) 
#(getting the ORDINARY generating function, accodring to the weight a^#edges of the set of
#labelled graphs on n vertices) and take
#the coefficient of a^(n-1) (the numbers we are interested in)
#That's the output
[seq(coeff(coeff(f,x,n)*n!,a,n-1),n=1..N)]:
end:


#ATreeSeq(N,r): The first N terms of the sequence eunerating the number of labelled "almost trees" on n vertices
#with the number of edges r more than that of a tree, i.e. n-1+r edges
#Using the edge-weighted egf for labelled connected graphs
#log(Sum((1+a)^binomial(n,2)*x^n/n!,i=0..infinity);
ATreeSeq:=proc(N,r) local f,x,a,n:

#Since we are only interested in the first N terms, we only need the terms up to x^(N) (the rest is irrelenat to us)
f:=log(add((1+a)^(n*(n-1)/2)*x^n/n!,n=0..N+1)):

#we take the Taylor exapansion up to x^N (the rest is garbage)
f:=taylor(f,x=0,N+1):

#We take the coefficient of x^n, multiply by n! (since it is an EXPONENTIAL generating function) 
#(getting the ORDINARY generating function, accodring to the weight a^#edges of the set of
#labelled graphs on n vertices) and take
#the coefficient of a^(n-1+r) (the numbers we are interested in)
#That's the output
[seq(coeff(coeff(f,x,n)*n!,a,n-1+r),n=1..N)]:
end:


#FunEqToSeq(PHI,z,N):
#Given a function PHI of z, and a positive integer N, outputs the first N terms, starting at n=1
#of the solutions of the functional equations
#f(x)=x*PHI(f(x)). Try:
#FunEqToSeq(exp(z),z,10);
FunEqToSeq:=proc(PHI,z,N) local f,i,j,x:

#we start with f=x, and keep applying the mapping f->x*PHI(f(x)), N times, and every time only
#save the terms up to and including x^N

f:=x:

for i from 1 to N do
f:=x*subs(z=f,PHI):
f:=taylor(f,x=0,N+1):
f:=add(coeff(f,x,j)*x^j,j=0..N):
od:

[seq(coeff(f,x,j),j=1..N)]:

end:


#LIF(PHI,z,N): The same as FunEqToSeq(PHI,z,N), but using the Lagrange Inversion Formula
#if f(x) is the solution of the functional equation
#f(x)=x*PHI(f(x)) then the coefficient of x^n in f(x) is 1/n times the coefficient of z^(n-1) in PHI(z)^n
LIF:=proc(PHI,z,N) local n:
[seq(coeff(taylor(PHI^n,z=0,n),z,n-1)/n,n=1..N)]:
end:


#TreeSeq1(N): Like TreeSeq(N), but using the fact that the number of ROOTED labeled trees if n! times
#the coefficient of the ogf T(x) that satisfies the functional equation
#T(x)=x*exp(T(x). Try:
#TreeSeq1(20);
TreeSeq1:=proc(N) local L,z,n:
L:=FunEqToSeq(exp(z),z,N):

[seq(L[n]*(n-1)!,n=1..N)]:
end:



#TreeSeqL(N,t): Like list of length N, whose  n-th entry is the (ordinary) generating function of
#ROTTED labelled trees on n vertices according to the weight t^#Leaves. It uses the fact that
#the exponential generating function satisfies the
#functional equation
#T(x,t)= t*x+ x*(exp(T(x,t)-1)
#Try
#TreeSeqL(20,t);
TreeSeqL:=proc(N,t) local L,z,n:
L:=FunEqToSeq(t-1+exp(z),z,N):

[seq(expand(L[n]*n!),n=1..N)]:
end:


#Added after class for the last homework problem for hw22
#AveAndMoms(f,x,K): given an ordinary generating function f of x, that is a weight-enumerator of some
#numerical attribute,and a positive integer K,
#outputs (i) The average  (ii) the standar-deviation, (iii) the 3rd through the K-th SCALED moments about the mean. Try
#AveAndMoms((1+x)^1000,x,6);
AveAndMoms:=proc(f,x,K) local f1,av,i,L,sd:
if subs(x=1,f)=0 then
 RETURN(FAIL):
fi:

f1:=f/subs(x=1,f):

av:=subs(x=1,diff(f1,x)):

f1:=f1/x^av:

f1:=x*diff(x*diff(f1,x),x):

sd:=sqrt(subs(x=1,f1)):

L:=[av,sd]:

for i from 3 to K do
f1:=x*diff(f1,x):
L:=[op(L),subs(x=1,f1)/sd^i]:
od:

evalf(L):

end: