######################################################################
##A435.txt: Save this file as  EnuTrees.txt                          #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read A435.txt<Enter>                               #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Rutgers University ,                  #
#zeilberg at math dot rutgers dot edu                                #
######################################################################
 
#Created: 
 
print(`Created:  June-July 2016`):
print(` This is A435.txt `):
print(`to automatically derive (rigorous!) expressions for the  first few moments of the random variable`):
print(`Total Height of labeled rooted trees`):
print(`it is inspired by the nice article`):
print(`The Enumeration of Rooted Trees by Total Height" by John Riordan and Neil Sloane`):
print(` J. Australian Math. Soc., 10 (1969), pp. 278-282.`):
print(`Available on line from: http://neilsloane.com/doc/riordan-enum-trees-by-height.pdf `):
print(` [viewed July 6, 2016] `):
print(``):
print(`It is the Maple package that accompanies the article `):
print(`Going Back to Neil Sloane's FIRST LOVE (OEIS Sequence A435): On the  Total Heights in Rooted Labeled Trees`):
print(`by Shalosh B. Ekhad and Doron Zeilberger`):
print(`published in the Personal Journal of Shalosh B. Ekhad and Doron Zeilberger`):
 print(`http://www.math.rutgers.edu/~zeilberg/pj.html  and arxiv.org .`):
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/  .`):

print(`---------------------------------------`):
 print(`For a list of the Supporting procedures type ezra1();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

print(`---------------------------------------`):
 print(`For a list of the MAIN procedures type ezra();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

with(combinat):

ezra1:=proc()

if args=NULL then
 print(` The supporting procedures are: Alpha, an, CheckJni, ChR, ExprFM, ExprM, Hafokh, FMomD, FMomL1, GP1, Jny, LimCOV, `):
 print(` Mom, MomAM, plotDist, RecJn, ScaledED, Zinn`):
 print(``):

else
ezra(args):
fi:

end:

ezra:=proc()

if args=NULL then
 print(`The main procedures are: aN, AlphaJ, Bitui, BituiNaked, ExprAM, ExprAMasy, FMom, FMomL, FMomF, Info, Jni, JNy, LimSM, LimSMa, Mamar,`):
 print(` MedianTH, MMM, ModeTH, SeqF, SeqFy, SeqFz, Summandk , Wn `):
 print(` `):


elif nops([args])=1 and op(1,[args])=Alpha then
print(`Alpha(f,x,N): Given a probability generating function`):
print(`f(x) (a polynomial, rational function and possibly beyond)`):
print(`returns a list, of length N, whose `):
print(`(i) First entry is the average`):
print(`(ii): Second entry is the variance`):
print(`for i=3...N, the i-th entry is the so-called alpha-coefficients`):
print(`that is the i-th moment about the mean divided by the`):
print(`variance to the power i/2 (For example, i=3 is the skewness`):
print(`and i=4 is the Kurtosis)`):
print(`If f(1) is not 1, than it first normalizes it by dividing`):
print(`by f(1) (if it is not 0) .`):
print(`For example, try:`):
print(`Alpha(((1+x)/2)^100,x,4);`):


elif nops([args])=1 and op(1,[args])=AlphaJ then
print(`AlphaJ(r,N): inputs a positive integer r, and a large positive integer N, outputs a list`):
print(`of length N, whose i-th item is a list of length r, consisting of (i) the average of the r.v."sum of heights" in a labeled rooted tree`):
print(`(ii) its variance, and the 3-rd through the r-th standardized moments (alpha coefficients). Try:`):
print(`AlphaJ(4,30);`):

elif nops([args])=1 and op(1,[args])=an then
print(`an(n): the  coefficient of x^n in R(x) times n!, where R(x) solves R(x)=x*exp(R(x)). Try: `):
print(`an(10); `):

elif nops([args])=1 and op(1,[args])=aN then
print(`aN(N): the first n coefficients of x^n in R(x) times n!, where R(x) solves R(x)=x*exp(R(x)). Try: `):
print(`aN(10); `):


elif nops([args])=1 and op(1,[args])=Bitui then
print(`Bitui(n,r,z): The expression in z and n such that`):
print(`its coefficient of z^n gives you  (d/dy)^r [J_n(y)] evaluated at y=1. Using Lagrange Inversion`):
print(`Try:`):
print(`Bitui(n,2,z);`):

elif nops([args])=1 and op(1,[args])=BituiNaked then
print(`BituiNaked(r,z): The expression in z such that  multiplied by (n-1)!*exp(n*z) gives you`):
print(`its coefficient of z^n gives you  (d/dy)^r [J_n(y)] evaluated at y=1. Using Lagrange Inversion`):
print(`Try:`):
print(`BituiNaked(2,z);`):

elif nops([args])=1 and op(1,[args])=CheckJni then
print(`CheckJni(K,N): Inputs  positive integers K and N, checks that the output of Jni(k,0) is what it should`):
print(`be for k from 1 to K, judging by the first N coefficients. Try:`):
print(`CheckJni(2,10); `):

elif nops([args])=1 and op(1,[args])=ChR then
print(`ChR(F,x,A,G): inputs an expression F the variable x, and  the symbol A , (denoting a function A of x), an expression G of A and x`):
print(`indicating that A(x) satisfies the first-order differential equation A'(x)=G(A(x),x)`):
print(`outputs the expression in A and x that d/dx (F(A(x),x)) equals to. Try:`):
print(`ChR(x*A^2,x, A,A/x/(1-A));`):


elif nops([args])=1 and op(1,[args])=ExprAM then
print(`ExprAM(r,n,W): the expression in n and W:=n!*add(n^k/k!,k=0..n-2) for`):
print(` the r-th moments about the  mean`):
print(` Try: `):
print( ` ExprAM(2,n,W); `):


elif nops([args])=1 and op(1,[args])=ExprAMasy then
print(`ExprAMasy(r,n): the leading term of the asymptotic expression in n for `):
print(` the r-th moments about the  mean of the r.v. "sum of heights of vertices" over all labeled rooted trees.`):
print(`Try: `):
print(`ExprAMasy(2,n);`):

elif nops([args])=1 and op(1,[args])=ExprFM then
print(`ExprFM(r,n,W): the expression in n and W:=n!*add(n^k/k!,k=0..n-2) for`):
print(` (d/dy)^r [J_n(y)] evaluated at y=1 (aka as factorial moment)`):
print(`Recall that by the Riordan-Sloane paper, W in asymptotic to n^n sqrt(Pi*n/2). `):
print(` Try: `):
print(` ExprFM(2,n,W); `):

elif nops([args])=1 and op(1,[args])=ExprM then
print(`ExprM(r,n,W): the expression in n and W:=n!*add(n^k/k!,k=0..n-2) for`):
print(` (y*d/dy)^r [J_n(y)] evaluated at y=1 (aka as straight moment)`):
print(`Recall that by the Riordan-Sloane paper, W in asymptotic to n^n sqrt(Pi*n/2). `):
print(` Try: `):
print(` ExprM(2,n,W); `):

elif nops([args])=1 and op(1,[args])=FMom then
print(`FMom(r,N):The first N terms, from n=1 to n=N,  of the sequence of r-th derivative of J_n(y) evaluated at y=1,`):
print(`done via the explicit expression given by Jni(n,0,R,x). Try:`):
print(`FMom(2,10);`):

elif nops([args])=1 and op(1,[args])=FMomF then
print(`FMomF(r,N): Like FMom(r,n) but using the formula giving by Summandk. Try:`):
print(`FMomF(2,10); `):

elif nops([args])=1 and op(1,[args])=FMomL then
print(`FMomL(r,N):The first N terms, from n=1 to n=N,  of the sequence of r-th derivative of J_n(y) evaluated at y=1,`):
print(`done via the explicit expression given by Jni(n,0,R,x) AND LAGRANGE INVERSION. Try:`):
print(`FMomL(2,10);`):

elif nops([args])=1 and op(1,[args])=FMomD then
print(`FMomD(r,N):The first N terms, from n=1 to n=N,  of the sequence of r-th derivative of J_n(y) evaluated at y=1,`):
print(`done directly. For checking purposes only. Try:`):
print(`FMomD(2,10);`):

elif nops([args])=1 and op(1,[args])=FMomL1 then
print(`FMomL1(r,n1): FMom(r,n)[n] using Bitui(n,r,z) that is gotten from Lagrange Inversion. Try:`):
print(`FMomL1(2,10):`):

elif nops([args])=1 and op(1,[args])=GP1 then
print(`GP1(L,k,d): guesses a polynomial of degree d fitting the list L, starting at k=1. Try:`):
print(`GP1([seq(i^2,i=1..10)],k,2);`):


elif nops([args])=1 and op(1,[args])=Hafokh then
print(`Hafokh(P,k,n): inputs a polynomial P in the variable k, of degree r, say, outputs the list of length r+1, consisting of `):
print(`polynomials in n, let's call it L`):
print(`such that`):
print(`P(k)=L[1]+L[2]*(n-k)+L[3]*(n-k)*(n-k-1)+...+L[r+1]*(n-k)*...*(n-k+r-1). Try:`):
print(`Hafokh(k^3,k,n);`):

elif nops([args])=1 and op(1,[args])=Info then
print(`Info(n,W,a,k): inputs symbols n and W, and a positive integer k, and outputs`):
print(`three lists concerning the random variable "sum of heights of vertices" defined on the`):
print(`set of n^(n-1) labeled rooted trees on n vertices. Where W stands for`):
print(` the expectation, given by n!/n^(n-1)*add(n^k/k!,k=0..n-2).`):
print(` The first is a list of length k whose first entry is W, and whose i-th entry is`):
print(`the explicit expression in terms of n, W, and Pi for the i-th moment about the mean`):
print(`The second list starts with [0,1] and for i from 2 to k is the limit, as n goes to infinity`):
print(`of the standardized moments, i.e. the moments of the limiting distibution.`):
print(`Try: `):
print(` Info(n,W,a,4); `):

elif nops([args])=1 and op(1,[args])=Jni then
print(`Jni(n,i,R,x): Inputs a non-negative integer n, a non-negative integer i, a symbol R, and a variable name x,`):
print(`outputs an expression, in R and x for the i-th derivative, with respect to x`):
print(`of J_n(x), in terms of R=J_0(x), where J_n(x) is the coefficient of z^n in the formal power series`):
print(`J(x,z) that satisfies the functional equation`):
print(`J(x,z)=x*exp(J(x+x*z,z)). Try:`):
print(`Jni(2,1,R,x);`):

elif nops([args])=1 and op(1,[args])=Jny then
print(`Jny(n,y): the coefficient of x^n in J(x) times n!, where J(x) solves J(x)=x*exp(J(x*y)). Try:`):
print(`Jny(5,y);`):

elif nops([args])=1 and op(1,[args])=JNy then
print(`JNy(N,y): the first n coefficients of x^n of J(x) times n!, where J(x) solves J(x)=x*exp(J(x*y)). Try:`):
print(`JNy(10,y); `):

elif nops([args])=1 and op(1,[args])=LimCOV then
print(`LimCOV(): the limit of the coefficient of variation of the random variable "total height" defined on labeled rooted trees on n vertices`):
print(`Try: `):
print(`LimCOV();`):

elif nops([args])=1 and op(1,[args])=LimSM then
print(`LimSM(r): the limit as n goes to infinity of the standardized`):
print(` r-th moments for  the r.v. "sum of heights of vertices" over all labeled rooted trees.`):
print(`Try: `):
print(` LimSM(4); `):

elif nops([args])=1 and op(1,[args])=LimSMa then
print(`LimSMa(r,a): the limit as n goes to infinity of the standardized`):
print(`r-th moments for  the r.v. "sum of heights of vertices" over all labeled rooted trees.`):
print(`in terms of the constant a=1/(10-3*Pi)`):
print(`Try: `):
print(`LimSMa(4,a);`):


elif nops([args])=1 and op(1,[args])=Mamar then
print(`Mamar(k): inputs a positive integer k and outputs`):
print(`an article about the moments, up to the k-th, of the random variable "sum of heights of vertices" defined on the`):
print(`set of n^(n-1) labeled rooted trees on n vertices. `):
print(`it gives explicit expressions, for the expectation, varaince, and the third through k-th moments (about the mean)`):
print(`it also gives the limit of the first k standardized moments as n goes to infiniy `):
print(`Try: `):
print(`Mamar(4);`):


elif nops([args])=1 and op(1,[args])=MedianTH then
print(`MedianTH(n): Out of the n^(n-1) rooted labeled trees, the median total height`):
print(`MedianTH(10);`):

elif nops([args])=1 and op(1,[args])=MMM then
print(`MMM(N): the sequences of (truncated) averages, modes, and median for the r.v. "total height" defined`):
print(`on rooted labeled trees, as well as the maximum possible total height. Try:`):
print(`MMM(20):`):

elif nops([args])=1 and op(1,[args])=ModeTH then
print(`ModeTH(n): Out of the n^(n-1) rooted labeled trees, what total height is the most likley, together with the number`):
print(`of trees that have that total height. Try:`):
print(`ModeTH(10);`):

elif nops([args])=1 and op(1,[args])=Mom then
print(`Mom(r,N):The first N terms, from n=1 to n=N,  of the sequence of of (y*d/dy)^r J_n(y) evaluated at y=1,`):
print(`done via Jni(r,0). Try:`):
print(`Mom(2,10);`):

elif nops([args])=1 and op(1,[args])=MomAM then
print(`MomAM(r,N):The first N terms, from n=1 to n=N,  of the sequence of of (y*d/dy)^r (J_n(y)/y^mean) evaluated at y=1,`):
print(`done via Jni(r,0). Try:`):
print(`MomAM(2,10);`):

elif nops([args])=1 and op(1,[args])=plotDist then
print(`plotDist(f,x,K1,K2): Given a prob. gen. function f(x) that has a `):
print(`Taylor series`):
print(`for a discrete r.v.`):
print(`plots its normalized version (X-mu)/sig between mu-K1*sig and`):
print(`mu+K2*sig`):
print(`For example, try:`):
print(`plotDist((1+x)^40,x,5,5);`):

elif nops([args])=1 and op(1,[args])=RecJn then
print(`RecJn(J,n,x):  Inputs a symbol J, a positive integer n, and a veriable name,x, outputs the`):
print(`expression  the coefficient of z^n, that we call J[n](x), in the solution of the functional equation`):
print(`J(x)=x*exp(J(x+x*z)) in terms of J[i](x) and its derivatives. The r-th derivative of J[i](x) is denoted by J[i,r](x)`):
print(`In particluar J[i,0](x) is J[i](x). Try:`):
print(`RecJn(J,2,x);`):

elif nops([args])=1 and op(1,[args])=ScaledED then
print(`ScaledED(f,x,K1,K2): Given a prob. gen. function f(x) that has a `):
print(`Taylor series,`):
print(`for a discrete r.v.`):
print(`outputs the  normalized version (X-mu)/sig between mu-K1*sig and`):
print(` mu+K2*sig `):
print(` For example, try: `):
print(`ScaledED((1+x)^40,x,3,3);`):

elif nops([args])=1 and op(1,[args])=SeqF then
print(`SeqF(F,x,N): inputs a function F representing a formal power series, `):
print(`phrased in terms of x, and outputs the first N coefficients of the solution to the Functional Equation`):
print(`R(x)=x*F(R(x)). Try:`):
print(` SeqF(exp(x),x,20); `):

elif nops([args])=1 and op(1,[args])=SeqFy then
print(`SeqFy(F,x,y,N): inputs a function F representing a formal power series, `):
print(`phrased in terms of x, and outputs the first N coefficients of the solution to the Functional Equation`):
print(`R(x)=x*F(R(x*y)). The output is a sequence of polynomials in y. Try:`):
print(` SeqFy(exp(x),x,y,10); `):

elif nops([args])=1 and op(1,[args])=SeqFz then
print(`SeqFz(F,x,z,k,N): inputs a function F representing a formal power series, `):
print(`phrased in terms of x, and outputs the first N coefficients, all truncated polynomials of z of degree k, of the solution to the Functional Equation`):
print(`R(x)=x*F(R(x*(1+z))). Try:`):
print(`SeqFz(exp(x),x,z,3,20);`):


elif nops([args])=1 and op(1,[args])=Summandk then
print(`Summandk(r,k): the polynomial in k, let's call it p(k), such that the`):
print(` (d/dy)^r [J_n(y)] evaluated at y=1 is given by the summation formula`):
print(`(n-1)!*add(p(k)*n^(n-k)/(n-k)!,k=0..n);`):
print(`Try:`):
print(`Summand(2,k);`):

elif nops([args])=1 and op(1,[args])=Wn then
print(`Wn(n): the sequence W_n defined in Eq. (3) of the Riordan-Sloane paper, namely OEIS sequence A1864, or n*A435`):
print(`n!*add(n^k/k!,k=0..n-2).  Try:`):
print(`Wn(10);`):
 
elif nops([args])=1 and op(1,[args])=Zinn then
print(`Zinn(resh): Zinn-Justin's method to estimate`):
print(`the C,mu, and theta such that`):
print(`resh[i] is appx. Const*mu^i*i^theta`):
print(`For example, try:`):
print(`Zinn([seq(5*i*2^i,i=1..30)]);`):

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

###start from Histabrut
with(plots):
#plotDist(f,x,K): Given a prob. gen. function f(x) that has a 
#Taylor series
#for a discrete r.v.
#plots its normalized version (X-mu)/sig between mu-K1*sig and
#mu+K2*sig
#For example, try:
#plotDist((1+x)^40,x,5,5);
plotDist:=proc(f,x,K1,K2) local mu,f1,lu,gu,sig,i,j1:
f1:=f/subs(x=1,f):
mu:=subs(x=1,diff(f1,x)):
gu:=f1/x^mu:
sig:=sqrt(subs(x=1,x*diff(x*diff(gu,x),x))):

lu:=taylor(f1,x=0,trunc(mu)+K2*trunc(sig)+10):

lu:=[seq([i,coeff(lu,x,i)],i=max(0,trunc(mu-K1*sig))..trunc(mu+K2*sig))]:

lu:=evalf([seq([(lu[j1][1]-mu)/sig,lu[j1][2]*sig],j1=1..nops(lu))]):

plot(lu,scaling=constrained):

end:

#ScaledED(f,x,K1,K2): Given a prob. gen. function f(x) that has a 
#Taylor series,
#for a discrete r.v.
#outputs the  normalized version (X-mu)/sig between mu-K1*sig and
#mu+K2*sig
#For example, try:
#ScaledED((1+x)^40,x,3,3);
ScaledED:=proc(f,x,K1,K2) local mu,f1,lu,gu,sig,i,j1:
f1:=f/subs(x=1,f):
mu:=subs(x=1,diff(f1,x)):
gu:=f1/x^mu:
sig:=sqrt(subs(x=1,x*diff(x*diff(gu,x),x))):

lu:=taylor(f1,x=0,trunc(mu)+K2*trunc(sig)+10):

lu:=[seq([i,coeff(lu,x,i)],i=max(0,trunc(mu-K1*sig))..trunc(mu+K2*sig))]:

lu:=evalf([seq([(lu[j1][1]-mu)/sig,lu[j1][2]*sig],j1=1..nops(lu))]):
lu:

end:

#Alpha(f,x,N): Given a probability generating function
#f(x) (a polynomial, rational function and possibly beyond)
#returns a list, of length N, whose 
#(i) First entry is the average
#(ii): Second entry is the variance
#for i=3...N, the i-th entry is the so-called alpha-coefficients
#that is the i-th moment about the mean divided by the
#variance to the power i/2 (For example, i=3 is the skewness
#and i=4 is the Kurtosis)
#If f(1) is not 1, than it first normalizes it by dividing
#by f(1) (if it is not 0) .
#For example, try:
#Alpha(((1+x)/2)^100,x,4);
Alpha:=proc(f,x,N) local gu,i:
 gu:=AveAndMoms(f,x,N):
 if gu=FAIL then
  RETURN(gu):
 fi:

if gu[2]=0 then
print(`The variance is 0`):
  RETURN(FAIL):
fi:

[gu[1],gu[2],seq(gu[i]/gu[2]^(i/2),i=3..N)]:

end:

#AveAndMoms(f,x,N): Given a probability generating function
#f(x) (a polynomial, rational function and possibly beyond)
#returns a list whose first entry is the average 
#(alias expectation, alias mean)
#followed by the variance, the third moment (about the mean) and
#so on, until the N-th moment (about the mean).
#If f(1) is not 1, than it first normalizes it by dividing
#by f(1) (if it is not 0) .
#For example, try:
#AveAndMoms(((1+x)/2)^100,x,4);

AveAndMoms:=proc(f,x,N) local mu,F,memu1,gu,i:
mu:=simplify(subs(x=1,f)):

if mu=0 then
print(f, `is neither a prob. generating function nor can it be made so`):
RETURN(FAIL):
fi:

F:=f/mu:


memu1:=simplify(subs(x=1,diff(F,x))):

gu:=[memu1]:

F:=F/x^memu1:

F:=x*diff(F,x):
for i from 2 to N do
 F:=x*diff(F,x):
 gu:=[op(gu),simplify(subs(x=1,F))]:
od:

gu:

end:

#From Zinn

sn:=proc(resh,n1):
-1/log(op(n1+1,resh)*op(n1-1,resh)/op(n1,resh)^2):
end:
 
#Zinn(resh): Zinn-Justin's method to estimate
#the C,mu, and theta such that
#resh[i] is appx. Const*mu^i*i^theta
#For example, try:
#Zinn([seq(5*i*2^i,i=1..30)]);
Zinn:=proc(resh)
local s1,s2,theta,mu,n1,i:
if nops({seq(sign(resh[i]),i=1..nops(resh))})<>1 then
 RETURN(FAIL):
fi:

n1:=nops(resh)-1:
s1:=sn(resh,n1):
s2:=sn(resh,n1-1):
theta:=evalf(2*(s1+s2)/(s1-s2)^2):
mu:=evalf(sqrt(op(n1+1,resh)/op(n1-1,resh))*exp(-(s1+s2)/((s1-s2)*s1))):
[theta,mu]:
end:

#End Zinn

###End from Histabrut


#SeqF(F,x,N): inputs a function F representing a formal power series, 
#phrased in terms of x, and outputs the first N coefficients of the solution to the Functional Equation
#R(x)=x*F(R(x)). Try:
#SeqF(exp(x),x,20);
SeqF:=proc(F,x,N) local F1,gu,i,i1:

if subs(x=0,F)=0 then
 RETURN(FAIL):
fi:

F1:=taylor(F,x=0,N+1):
F1:=add(coeff(F1,x,i)*x^i,i=0..N):

gu:=x*subs(x=0,F):
for i from 1 to N do
gu:=x*subs(x=gu,F1):
gu:=taylor(gu,x=0,N+1):
gu:=add(coeff(gu,x,i1)*x^i1,i1=0..N):
od:

[seq(coeff(gu,x,i1),i1=1..N)]:
end:

#SeqFy(F,x,y,N): inputs a function F representing a formal power series, 
#phrased in terms of x, and outputs the first N coefficients, all polynomials of y, of the solution to the Functional Equation
#R(x)=x*F(R(x*y)). Try:
#SeqF(exp(x),x,y,20);
SeqFy:=proc(F,x,y,N) local F1,gu,i,i1:

if subs(x=0,F)=0 then
 RETURN(FAIL):
fi:

F1:=taylor(F,x=0,N+1):
F1:=add(coeff(F1,x,i)*x^i,i=0..N):

gu:=x*subs(x=0,F):
for i from 1 to N do
gu:=x*subs(x=subs(x=x*y,gu),F1):
gu:=taylor(gu,x=0,N+1):
gu:=add(coeff(gu,x,i1)*x^i1,i1=0..N):
od:

[seq(coeff(expand(gu),x,i1),i1=1..N)]:
end:


#an(n): the coefficient of x^n in R(x) times n! where R(x) solves R(x)=x*exp(R(x))
an:=proc(n) local k:
option remember:
if n=1 then
 1:
else
add(k*an(k)*an(n-k),k=1..n-1)/(n-1):
fi:
end:

#aN(N): the first n coefficients of x^n in R(x) times n! where R(x) solves R(x)=x*exp(R(x))
aN:=proc(N) local n:
[seq(n!*an(n),n=1..N)]:
end:


#Jny(n,y): the coefficient of x^n of J(x) where J(x) solves J(x)=x*exp(J(x*y))
Jny:=proc(n,y) local k:
option remember:
if n=1 then
 1:
else
expand(add(k*Jny(k,y)*Jny(n-k,y)*y^k,k=1..n-1)/(n-1)):
fi:
end:

#JNy(N,y): the first n coefficients of x^n of J(x) where J(x) solves J(x)=x*exp(J(x*y))
JNy:=proc(N,y) local n:
[seq(expand(n!*Jny(n,y)),n=1..N)]:
end:


#SeqFz(F,x,z,k,N): inputs a function F representing a formal power series, 
#phrased in terms of x, and outputs the first N coefficients, all truncated polynomials of z of degree k, of the solution to the Functional Equation
#R(x)=x*F(R(x*(1+z))). Try:
#SeqFz(exp(x),x,z,3,20);
SeqFz:=proc(F,x,z,k,N) local F1,gu,i,i1:

if subs(x=0,F)=0 then
 RETURN(FAIL):
fi:

F1:=taylor(F,x=0,N+1):
F1:=add(coeff(F1,x,i)*x^i,i=0..N):

gu:=x*subs(x=0,F):
for i from 1 to N do
gu:=x*subs(x=subs(x=x*(1+z),gu),F1):
gu:=taylor(gu,z=0,k+1):
gu:=add(coeff(gu,z,i1)*z^i1,i1=0..k):
gu:=taylor(gu,x=0,N+1):
gu:=add(coeff(gu,x,i1)*x^i1,i1=0..N):
od:

[seq(coeff(expand(gu),x,i1),i1=1..N)]:
end:

#RecJn(J,n,x):  Inputs a symbol J, a positive integer n, and a veriable name,x, outputs the
#expression  the coefficient of z^n, that we call J[n](x), in the solution of the functional equation
#J(x)=x*exp(J(x+x*z)) in terms of J[i](x) and its derivatives. The r-th derivative of J[i](x) is denoted by J[i,r](x)
#In particluar J[i,0](x) is J[i](x). Try:
#RecJn(J,2,x);
RecJn:=proc(J,n,x) local i,gu,m,z,c:
gu:=add(z^m*add(x^(m-i)/(m-i)!*J[i,m-i],i=0..m),m=0..n):
gu:=x*exp(gu):
gu:=taylor(gu,z=0,n+1):
gu:=normal(coeff(gu,z,n)):
gu:=normal(subs(exp(J[0,0])=J[0,0]/x,gu)):
c:=coeff(gu,J[n,0],1):
gu:=expand(gu-c*J[n,0]):
gu/(1-c):
end:

#ChR(F,x,A,G): inputs an expression F the variable x, and  the symbol A , (denoting a function A of x), an expression G of A and x
#indicating that A(x) satisfies the first-order differential equation A'(x)=G(A(x),x)
#outputs the expression in A and x that d/dx (F(A(x),x)) equals to. Try
#ChR(x*A^2,x, A,A/x/(1-A));
ChR:=proc(F,x,A,G) local gu,f:
gu:=subs(A=f(x),F):
gu:=diff(gu,x):
gu:=subs(diff(f(x),x)=G,gu):
gu:=subs(f(x)=A,gu):
normal(gu):
end:



#Jni(n,i,R,x): Inputs a non-negative integer n, a non-negative integer i, a symbol R, and a variable name x,
#outputs an expression, in R and x for the i-th derivative, with respect to x
#of J_n(x), in terms of R=J_0(x), where J_n(x) is the coefficient of z^n in the formal power series
#J(x,z) that satisfies the functional equation
#J(x,z)=x*exp(J(x+x*z,z)). Try
#Jni(2,1,R,x);
Jni:=proc(n,i,R,x) local gu,J,i1,n1:
option remember:

if not (type(n,integer) and n>=0 and type(i,integer) and i>=0 and type(R,symbol) and type(x,symbol)) then
 RETURN(FAIL):
fi:

if n=0 and i=0 then
 RETURN(R):
fi:

if i>0 then
 gu:=Jni(n,i-1,R,x):
 gu:=ChR(gu,x,R,R/x/(1-R)):
 RETURN(gu):
fi:

gu:=RecJn(J,n,x):

for n1 from 0 to n-1 do
for i1 from 0 to n+1 do
gu:=normal(subs(J[n1,i1]=Jni(n1,i1,R,x),gu)):
od:
od:

gu:

end:



#CheckJni(K,N): Inputs  positive integers K and N, checks that the output of Jni(k,0) is what it should
#be for k from 1 to K, judging by the first N coefficients. Try:
#CheckJni(2,10);
CheckJni:=proc(N,K) local R,x,gu,i1,R1,lu,R0,z,k:
gu:=SeqFz(exp(x),x,z,K,N):
R0:=add(coeff(gu[i1],z,0)*x^i1,i1=1..N):

for k from 1 to K do
 R1:=add(coeff(gu[i1],z,k)*x^i1,i1=1..N):
 lu:=Jni(k,0,R,x):
 lu:=subs(R=R0,lu)-R1:
 lu:=taylor(lu,x=0,N+1):
 if {seq(coeff(lu,x,i1),i1=0..N)}<{0} then
   print(`For k=`, k, `something is wrong `):
  RETURN(false):
 fi:
od:

true:

end:




#FMomD(r,N):The first N terms, from n=1 to n=N,  of the sequence of r-th derivative of J_n(y) evaluated at y=1,
#done directly. For checking purposes only. Try:
#FMomD(2,10);
FMomD:=proc(r,N) local gu,y,i:
gu:=JNy(N,y):
if r>0 then
gu:=[seq(diff(gu[i],y$r),i=1..N)]:
fi:

gu:=subs(y=1,gu):
end:

#FMom(r,N):The first N terms, from n=1 to n=N,  of the sequence of r-th derivative of J_n(y) evaluated at y=1,
#done via Jni(r,0). Try:
#FMom(2,10);
FMom:=proc(r,N) local gu,x,R0,i,R:
option remember:
gu:=Jni(r,0,R,x):
R0:=add(i^(i-1)/i!*x^i,i=1..N):
gu:=r!*subs(R=R0,gu):
gu:=taylor(gu,x=0,N+1):

[seq(i!*coeff(gu,x,i),i=1..N)]:

end:





#Mom(r,N):The first N terms, from n=1 to n=N,  of the sequence of of (y*d/dy)^r J_n(y) evaluated at y=1,
#done via Jni(r,0). Try:
#Mom(2,10);
Mom:=proc(r,N) local gu,i,j:
option remember:

gu[0]:=[seq(i^(i-1),i=1..N)]:

for i from 1 to r do
 gu[i]:=FMomF(i,N):
od:

[seq(add(stirling2(r,i)*gu[i][j],i=0..r),j=1..N)]:

end:

#MomAM(r,N):The first N terms, from n=1 to n=N,  of the sequence of of (y*d/dy)^r (J_n(y)/y^mean) evaluated at y=1,
#done via Jni(r,0). Try:
#MomAM(2,10);
MomAM:=proc(r,N) local gu0,gu,i,j:
option remember:

gu0:=FMom(0,N):

for i from 0 to r do
 gu[i]:=Mom(i,N):
 gu[i]:=[seq(gu[i][j]/gu0[j],j=1..N)]:
od:

[seq((-gu[1][j])^r+add(binomial(r,i)*gu[i][j]*(-gu[1][j])^(r-i),i=1..r),j=1..N)]:

end:

#AlphaJ(r,N): inputs a positive integer r, and a large positive integer N, outputs a list
#of length N, whose i-th item is a list of length r, consisting of (i) the average of the r.v."sum of heights" in a labeled rooted tree
#(ii) its variance, and the 3-rd through the r-th standardized moments (alpha coefficients). Try:
#Alpha(4,30);
AlphaJ:=proc(r,N) local gu0,gu,i,j:

gu0:=Mom(0,N):

gu[1]:=Mom(1,N):

gu[1]:=[seq(gu[1][j]/gu0[j],j=1..N)]:

for i from 2 to r do
gu[i]:=MomAM(i,N):
od:

for i from 3 to r do
 gu[i]:=[FAIL,FAIL,seq(gu[i][j]/gu[2][j]^(i/2),j=3..N)]
od:

[seq([seq(gu[i][j],i=1..r)],j=1..N)]:
end:





#Bitui(n,r,z): The expression in z and n such that
#its coefficient of z^n gives you  (d/dy)^r [J_n(y)] evaluated at y=1. Using Lagrange Inversion
#Try:
#Bitui(n,2,z);
Bitui:=proc(n,r,z) local gu,x:
gu:=normal(diff(Jni(r,0,z,x),z)):

(n-1)!*z*exp(n*z)*gu:
end:


#FMomL1(r,n1): FMom(r,n)[n] using Bitui(n,r,z). Try:
#FMomL1(2,10):
FMomL1:=proc(r,n1) local gu,n,z:
gu:=Bitui(n,r,z):
gu:=subs(n=n1,gu):
eval(r!*coeff(taylor(gu,z=0,n1+1),z,n1)):
end:


#FMomL(r,N): Like FMom(r,n) using Bitui(n,r,z). Try:
#FMomL(2,10):
FMomL:=proc(r,N) local n1:
[seq(FMomL1(r,n1),n1=1..N)]:
end:



#BituiNaked(r,z): The expression in z such that  multiplied by (n-1)!*exp(n*z) gives you
#its coefficient of z^n gives you  (d/dy)^r [J_n(y)] evaluated at y=1. Using Lagrange Inversion
#Try:
#BituiNaked(2,z);
BituiNaked:=proc(r,z) local gu,x:
gu:=normal(diff(Jni(r,0,z,x),z)):

z*gu:

end:

#GP1(L,k,d): guesses a polynomial of degree d fitting the list L, starting at k=1. Try:
#GP1([seq(i^2,i=1..10)],k,2);
GP1:=proc(L,k,d) local eq,var,a,i,P:
if nops(L)-d<3 then
  RETURN(FAIL):
fi:

var:={seq(a[i],i=0..d)}:
P:=add(a[i]*k^i,i=0..d):
eq:={seq(subs(k=i,P)-L[i],i=1..d+2)}:
var:=solve(eq,var):
if var=NULL then
 RETURN(FAIL):
fi:
P:=subs(var,P):
if {seq(subs(k=i,P)-L[i],i=1..nops(L))}<>{0} then
  RETURN(FAIL):
fi:
P:
end:



#Summandk(r,k): the polynomial in k, let's call it p(k), such that the
# (d/dy)^r [J_n(y)] evaluated at y=1 is given by the summation formula
#(n-1)!*add(p(k)*n^(n-k)/(n-k)!,k=0..n);
#Try:
#Summand(2,k);
Summandk:=proc(r,k) local gu,z,d,i:
gu:=factor(BituiNaked(r,z)):

d:=degree(denom(gu),z)-1:

gu:=taylor(gu,z=0,d+11):
gu:=[seq(coeff(gu,z,i),i=1..d+10)]:
gu:=GP1(gu,k,d):
r!*gu:
end:



#FMomF(r,N): Like FMom(r,n) but using the formula giving by Summandk. Try:
#FMomF(2,10);
FMomF:=proc(r,N) local P,k,k1,n1:

if r=0 then
 RETURN(FMom(0,N)):
fi:
P:=Summandk(r,k):

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

[seq((n1-1)!*add(subs(k=k1,P)*n1^(n1-k1)/(n1-k1)!,k1=0..n1),n1=1..N)]:

end:


#Hafokh(P,k,n): inputs a polynomial P in the variable k, of degree r, say, outputs the list of length r+1, consisting of 
#polynomials in n, let's call it L
#such that
#P(k)=L[1]+L[2]*(n-k)+L[3]*(n-k)*(n-k-1)+...+L[r+1]*(n-k)*...*(n-k+r-1). Try:
#Hafokh(k^3,k,n);
Hafokh:=proc(P,k,n) local r,i,eq,var,a,gu,j,L:
r:=degree(P,k):
var:={seq(a[i],i=0..r)}:

gu:=expand(add(a[i]*mul(n-k-j,j=0..i-1),i=0..r)-P):
eq:={seq(coeff(gu,k,i),i=0..r)}:

var:=solve(eq,var):

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

L:=[seq(subs(var,a[i]),i=0..r)]:

if expand(add(L[i+1]*mul(n-k-j,j=0..i-1),i=0..r)-P)<>0 then
 FAIL:
else
 L:
fi:

end:



#ExprFM(r,n,W): the expression in n and W:=n!*add(n^k/k!,k=0..n-2) for
# (d/dy)^r [J_n(y)] evaluated at y=1 (aka as factorial moment)
#Try:
#ExprFM(2,n,W);
ExprFM:=proc(r,n,W) local lu,gu,k,i,j:
option remember:
if r=0 then
 RETURN(1):
fi:

lu:=Hafokh(Summandk(r,k),k,n):

gu:=0:

for i from 0 to nops(lu)-1 do
 gu:=gu+lu[i+1]*n^i*(W-add(n^(n-j)*simplify((n-1)!/(n-j)!),j=0..i-1)):
od:

gu:=expand(simplify(subs(W=W/n+2*n^(n-1),gu ))):
gu:=subs(W=W*n^(n-1),gu):
gu:=simplify(gu/n^(n-1)):
expand(gu):

end:


#ExprM(r,n,W): the expression in n and W:=n!*add(n^k/k!,k=0..n-2) for
# (y*d/dy)^r [J_n(y)] evaluated at y=1 (aka as straight moments)
#Try:
#ExprM(2,n,W);
ExprM:=proc(r,n,W) local gu,i:
option remember:

for i from 0 to r do
 gu[i]:=ExprFM(i,n,W):
od:

add(stirling2(r,i)*gu[i],i=0..r):

end:

#ExprAM(r,n,W): the expression in n and W:=n!*add(n^k/k!,k=0..n-2) for
# the r-th moments about the  mean
#Try:
#ExprAM(2,n,W);
ExprAM:=proc(r,n,W) local gu,i,i1:
option remember:
for i from 1 to r do
 gu[i]:=ExprM(i,n,W):
od:
gu:=expand((-gu[1])^r+add(binomial(r,i1)*gu[r-i1]*(-gu[1])^i1,i1=0..r-1)):
collect(gu,W):
end:




#ExprAMasy(r,n): the leading term of the asymptotic expression in n for 
# the r-th moments about the  mean of the r.v. "sum of heights of vertices" over all labeled rooted trees.
#Try:
#ExprAMasy(2,n,W);
ExprAMasy:=proc(r,n) local gu,m,W,d:
option remember:
gu:=ExprAM(r,n,W):
gu:=expand(subs({n=m^2,W=sqrt(Pi/2)*m^3},gu));
d:=degree(gu,m):
coeff(gu,m,d)*n^(d/2):
end:


#LimSM(r): the limit as n goes to infinity of the standardized
# r-th moments for  the r.v. "sum of heights of vertices" over all labeled rooted trees.
#Try:
#LimSM(4);
LimSM:=proc(r) local n:
option remember:
if n mod 2=0 then
normal(ExprAMasy(r,n)/ExprAMasy(2,n)^(r/2)):
else
sqrt(normal(ExprAMasy(r,n)^2/ExprAMasy(2,n)^r)):
fi:
end:

#LimSMa(r,a): the limit as n goes to infinity of the standardized
# r-th moments for  the r.v. "sum of heights of vertices" over all labeled rooted trees.
#in terms of the constant a=1/(10-3*Pi)
#Try:
#LimSMa(4,a);
LimSMa:=proc(r,a) local gu:
option remember:
gu:=LimSM(r):
normal(subs(Pi=(10*a-1)/3/a,gu)):


end:

#Info(n,W,a,k): inputs symbols n and W, a, and a positive integer k, and outputs
#three lists concerning the random variable "sum of heights of vertices" defined on the
#set of n^(n-1) labeled rooted trees on n vertices. Where W stands for
# the expectation, given by n!/n^(n-1)*add(n^k/k!,k=0..n-2).
# The first is a list of length k whose first entry is W, and whose i-th entry is
#the explicit expression in terms of n, W, and Pi for the i-th moment about the mean
#The second list starts with [0,1] and for i from 2 to k is the limit, as n goes to infinity
#of the standardized moments, i.e. the moments of the limiting distibution.
#The third is the rephrasing the second in terms of the constant a=1/(10-3*Pi), since
#it turns out that these standardized moments are polynomials in a (or square-roots of them for odd moments)
#Try:
#Info(n,W,a,4);
Info:=proc(n,W,a,k) local i:
option remember:
[[W,seq(ExprAM(i,n,W),i=2..k)],
[0,1,seq(LimSM(i),i=3..k)],
[0,1,seq(LimSMa(i,a),i=3..k)]]:
end:



#Mamar(k): inputs a positive integer k and outputs
#an article about the moments, up to the k-th, of the random variable "sum of heights of vertices" defined on the
#set of n^(n-1) labeled rooted trees on n vertices. 
# the expectation, given by n!/n^(n-1)*add(n^k/k!,k=0..n-2).
# The first is a list of length k whose first entry is W, and whose i-th entry is
#the explicit expression in terms of n, W, and Pi for the i-th moment about the mean
#The second list starts with [0,1] and for i from 2 to k is the limit, as n goes to infinity
#of the standardized moments, i.e. the moments of the limiting distibution.
#The third is the rephrasing the second in terms of the constant a=1/(10-3*Pi), since
#it turns out that these standardized moments are polynomials in a (or square-roots of them for odd moments)
#Try:
#Mamar(4);
Mamar:=proc(k) local n,V,W,a,x,R,J,y,i,gu,t0:
t0:=time():
print(``):
print(`Explicit (and Asymptotic) Expressions for the Expectation, Variance, and all Moments (about the mean) up to the`, k, `-th `):
print(`of the random variable: Total Height defined on the set of Rooted Labeled Trees with n vertices `):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`Arthur Cayley famously proved that the number of labeled rooted trees on n vertices, let's call it r(n),  equals`, n^(n-1), ` . `):
print(`One of the many ways of proving it is by noting that the exponential generating function, let's call it`, R(x), `defined by` ):
print(``):
print(R(x)=Sum(r(n)*x^n/n!,n=0..infinity)):
print(``):
print(`satisfies the Functional Equation, that follows from standard generatingfunctionlogy `):
print(``):
print(R(x)=x*exp(R(x)) ):
print(``):
print(`and it follows from the famous Lagrange Inversion Formula (e.g. http://www.math.rutgers.edu/~zeilberg/mamarim/mamarimPDF/lag.pdf)`):
print(``):
print(`that `, r(n)/n!, `equals  1/n times the coefficient of `, x^(n-1), `in the formal power series`, (exp(x))^n):
print(`hence r(n) equals (n-1)! times n^(n-1)/(n-1)!= n^(n-1) `, ` . `):
print(``):
print(`In a delightful article (and very historically significant, since that's how the OEIS started!)`):
print(` [Available on-line from http://neilsloane.com/doc/riordan-enum-trees-by-height.pdf], entitled `):
print(` "The Enumeration of Rooted Trees by Total Height"`):
print(` J. Australian Math. Soc., 10 (1969), pp. 278-282,`):
print(`  John Riordan and Neil Sloane considered the polynomials`, J[n](y), `defined as the weight-enumerators of the set of rooted labeled trees `):
print(`on n vertices, according to the weight, "total height", where the height of a vertex is its distance to the root`):
print(`and the total height of a rooted labeled tree is the sum of the heights (distance to the root) of the vertices`):
print(`Let J(x,y) be the exponential generating function of the sequence J[n](y)`):
print(``):
print(J(x,y)=Sum(J[n](y)*x^n/n!,n=0..infinity), `. `):
print(``):
print(`The same generatingfunctionology argument that lead to the functional equation for R(x), leads to the functional equation`):
print(``):
print(J(x,y)=x*exp(J(x*y,y))):
print(``):
print(`Alas, it does not seem to be possible to get a closed-form expression for the polynomials  J[n](y) (for symbolic n) . `):
print(``):
print(`For the record, here are the first 12 members of the sequence of polynomials`, J[n](y) ):
gu:=JNy(12,y):
print(``):
print(gu):
print(``):
print(`and in Maple notation`):
print(``):
lprint(gu):
print(``):

print(`Nevertheless, it is still interesting to ivestigate the statistical properties of the random variable`):
print(`total height, defined on the "sample space" of rooted labeled trees on n vertices.`):
print(`Note that`, J[n](1)=n^(n-1) , `and the sum of all total heights (that would yield the expectation upon dividing by`, n^(n-1) ):
print(`is the FIRST derivative of`, J[n](y), `with respect to y, evaluated at y=1. `):
print(`This is sequence https://oeis.org/A001864 in the OEIS,`):
print(`whose first 14 terms are: `):
print(0, 2, 24, 312, 4720, 82800, 1662024, 37665152, 952401888, 26602156800, 813815035000, 27069937855488, 972940216546896, 37581134047987712):
print(``):
print(`Divided by n (to get the mean height among all vertices) it is https://oeis.org/A000435 `):
print(``):
print(` 0, 1, 8, 78, 944, 13800, 237432, 4708144, 105822432, 2660215680, 73983185000, 2255828154624, 74841555118992, 2684366717713408, 103512489775594200`, `, `):
print(``):
print(`that according to the comment there (written by Neil Sloane) it is "The sequence that started it all" `):
print(``):
print(`John Riordan and Neil Sloane found an explicit formula for that sequence`):
print(``):
print(V[n]= n!*Sum(n^i/i!,i=0..n-2) ):
print(``):
print(`and then used G.N. Watson's elaboration of a formula of Ramanujan to deduce that the average total height, that we will call `, W[n]):
print(W[n]= (n!/n^(n-1))  *Sum(n^i/i!,i=0..n-2)):
print(``):
print(`i.e. `, W[n]=V[n]/n^(n-1), ` is asymptotic to`):
print(``):
print(sqrt(Pi/2)*n^(3/2)):
print(``):
print(`Hence the asymptotic mean height of random rooted labled tree is`, sqrt(Pi/2)*sqrt(n)):
print(`and in Maple notation`):
lprint(sqrt(Pi/2)*sqrt(n)):
print(`The purpose of this article is to state (rigorously proved!) explicit expressions for the variance, and the moments (about the mean)`):
print(`up to the`, k, `-th, in terms of n and W[n] `):
print(`as well as their asymptotics (implied by the Watson-Ramanujan formula)`):
print(`and the limits of the standardized moments, that would hopefully lead, one day, to an explicit evaluation of the probability`):
print(`density function of the scaled-limiting distribution. `):

print(``):
print(`---------------------------------------------------------------------------------------------`):
print(``):

print(`Theorem 1 (Riordan-Sloane): The expectation of the "total height" in the set of rooted labeled trees is`):
print(``):
print(W[n]):
print(``):
print(`whose asympotitc expression is`, sqrt(Pi/2)*n^(3/2) ):
print(``):
print(`and in Maple notation`):
lprint(sqrt(Pi/2)*n^(3/2)):
print(``):
print(`All the remaining`, k-1 , `theorem are new, as far as I know. In the statements below`):
print(a=1/(10-3*Pi)):
print(``):
print(`and in Maple notation`):
print(``):
lprint(a=1/(10-3*Pi)):

print(``):


if k>=2 then
print(``):
print(`---------------------------------------------------------------------------------------------`):
print(``):

print(`Theorem 2: The variance of the random variable "total height" on the set of rooted labeled trees on n vertices is given explicitly by`):
gu:= subs(W=W[n],ExprAM(2,n,W)):
print(``):
print(gu):
print(``):
print(`and in Maple notation`):
print(``):
lprint(gu):
print(``):

gu:=ExprAMasy(2,n):
print(`its asymptotic expression is`):
print(``):
print(gu):
print(``):
print(`and in Maple notation`):
print(``):
lprint(gu):
print(``):
gu:=normal(gu/n^3):
gu:=sqrt(gu)/sqrt(Pi/2):

print(`Hence the limit of the coefficient of variation, as n goes to infinity is`, gu):
print(``):
print(`and in Maple notation`):
print(``):
lprint(gu):
print(``):
print(`that is approximately equal to`, Digits, `significant digits to`):
print(``):
print(evalf(gu)):
print(``):
print(`so the limiting coefficient of variation is`, evalf(100*gu), `percents, and there is NO concentration about the mean.`):

fi:




if k>=3 then
print(``):
print(`---------------------------------------------------------------------------------------------`):
print(``):

print(`Theorem 3: The third moment (about the  mean) of the `):
print(`random variable "total height" on the set of rooted labeled trees on n vertices is given explicitly by`):
gu:= subs(W=W[n],ExprAM(3,n,W)):
print(``):
print(gu):
print(``):
print(`and in Maple notation`):
print(``):
lprint(gu):
print(``):

gu:=ExprAMasy(3,n):
print(`its asymptotic expression is`):
print(``):
print(gu):
print(``):
print(`and in Maple notation`):
print(``):
lprint(gu):
print(``):

print(`in floating point`):
print(evalf(gu)):
print(``):
print(`and in Maple notation`):
print(``):
lprint(evalf(gu)):
print(``):

gu:=LimSM(3):

print(`Hence the limit of the SKEWNESS, as n goes to infinity is`, gu):
print(``):
print(`that in Maple notation is`):
print(``):
print(lprint(gu)):
print(``):
print(`this is approximately to`, Digits, `digits `):
print(``):
print(lprint(evalf(gu))):
print(``):
gu:=LimSMa(3,a):
print(`and in terms of  a=1/(10-3*Pi) this equals `,gu):
print(``):
print(`that in Maple notation is`):
lprint(gu):
print(``):
print(`hence the limiting distribution is positively skewed, `):
fi:


if k>=4 then
print(``):
print(`---------------------------------------------------------------------------------------------`):
print(``):

print(`Theorem 4: The fourth moment (about the  mean) of the `):
print(`random variable "total height" on the set of rooted labeled trees on n vertices is given explicitly by`):
gu:= subs(W=W[n],ExprAM(4,n,W)):
print(``):
print(gu):
print(``):
print(`and in Maple notation`):
print(``):
lprint(gu):
print(``):

gu:=ExprAMasy(4,n):
print(`its asymptotic expression is`):
print(``):
print(gu):
print(``):
print(`and in Maple notation`):
print(``):
lprint(gu):
print(``):

print(`in floating point`):
print(evalf(gu)):
print(``):
print(`and in Maple notation`):
print(``):
lprint(evalf(gu)):
print(``):

gu:=LimSM(4):

print(`Hence the limit of the KURTOSIS, as n goes to infinity is`, gu):
print(``):
print(`that in Maple notation is`):
print(``):
print(lprint(gu)):
print(``):
print(`this is approximately to`, Digits, `digits `):
print(``):
print(lprint(evalf(gu))):
print(``):
gu:=LimSMa(4,a):
print(`and in terms of  a=1/(10-3*Pi) this equals `,gu):
print(``):
print(`that in Maple notation is`):
lprint(gu):
print(``):
print(`hence the limiting distribution is  LEPTOKURTIC. `):

fi:




for i from 5 to k do
print(``):
print(`---------------------------------------------------------------------------------------------`):
print(``):

print(`Theorem Number`, i, ` The`, i, `-th moment (about the mean) of the random variable "total height" on the set of rooted labeled trees`):
print(` on n vertices is given explicitly by`):
gu:= subs(W=W[n],ExprAM(i,n,W)):
print(``):
print(gu):
print(``):
print(`and in Maple notation`):
print(``):
lprint(gu):
print(``):

gu:=ExprAMasy(i,n):
print(`its asymptotic expression is`):
print(``):
print(gu):
print(``):
print(`and in Maple notation`):
print(``):
lprint(gu):
print(``):

print(`in floating point`):
print(evalf(gu)):
print(``):
print(`and in Maple notation`):
print(``):
lprint(evalf(gu)):
print(``):

gu:=LimSM(i):

print(`Hence the limit of the scaled`, i, `-th momemnt , as n goes to infinity is`, gu):
print(``):
print(`that in Maple notation is`):
print(``):
print(lprint(gu)):
print(``):
print(`this is approximately to`, Digits, `digits `):
print(``):
print(lprint(evalf(gu))):
print(``):
gu:=LimSMa(i,a):
print(`and in terms of  a=1/(10-3*Pi) this equals `,gu):
print(``):
print(`that in Maple notation is`):
lprint(gu):
print(``):
od:

print(`To sum up, the limits of the scaled momentes (aka alpha coefficiets, from the 3-rd to the`, k, `-th are`):
gu:=[seq(LimSM(i),i=3..k)]:
print(``):
print(gu):
print(``):
print(`and in Maple notation:`):
print(``):
lprint(gu):
print(``):
print(`and in floating point:`):
print(``):
lprint(evalf(gu)):
print(``):
gu:=[seq(LimSMa(i,a),i=3..k)]:
print(``):
print(`and in terms of a=1/(10-3*Pi) `):
print(``):
print(gu):
print(``):
print(`and in Maple notation:`):
print(``):
lprint(gu):
print(``):
print(`--------------------------------------------------------------`):
print(``):

print(`This ends this article that took`, time()-t0, `seconds to generate. `):

end:



#LimCOV(): the limit of the coefficient of variation
LimCOV:=proc() local gu,n:
gu:=ExprAMasy(2,n):
gu:=normal(gu/n^3):
gu:=sqrt(gu)/sqrt(Pi/2):
end:


#Wn(n): the sequence W_n defined in Eq. (3) of the Riordan-Sloane paper, namely
#n!*add(n^k/k!,k=0..n-2).  Try:
#Wn(10);
Wn:=proc(n) local k:
n!*add(n^k/k!,k=0..n-2):
end:

#WnAsy(n): the asymptotic formula for Wn(n)/n^(n). Try
#WnAsy(1000):
WnAsy:=proc(n)
sqrt(Pi/2)*sqrt(n)-4/3+sqrt(Pi/2)/12/sqrt(n)-4/135/n:
end:


#ModeTH(n): Out of the n^(n-1) rooted labeled trees, what total height is the most likley, together with the number
#of trees that have that total height. Try:
#ModeTH(10);
ModeTH:=proc(n) local gu,y,i,si,aluf:
gu:=n!*Jny(n,y):
aluf:=1:
si:=coeff(gu,y,1):

for i from 2 to degree(gu,y) do
if coeff(gu,y,i)>si then
 si:=coeff(gu,y,i):
 aluf:=i:
fi:
od:

[aluf,si]:

end:

#MedianTH(n): Out of the n^(n-1) rooted labeled trees, the median total height
#MedianTH(10);
MedianTH:=proc(n) local gu,gu0,i,su,y:

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

gu:=n!*Jny(n,y):


gu0:=subs(y=1,gu)/2:
su:=0:

for i from 1 to degree(gu,y) do
su:=su+coeff(gu,y,i):
if su>gu0 then
 RETURN(i):
fi:
od:
FAIL:
end:



#MMM(N): the sequences of (truncated) averages, modes, and median for the r.v. "total height" defined
#on rooted labeled trees, as well as the maximum possible total height. Try:
#MMM(20):
MMM:=proc(N) local i,gu1,gu2,gu3,gu2a,gu2b:

gu1:=[seq(trunc(Wn(i)/i^(i-1)),i=1..N)]:

gu2:=[seq(ModeTH(i),i=1..N)]:

gu2a:=[seq(gu2[i][1],i=1..nops(gu2))]:
gu2b:=[seq(gu2[i][2],i=1..nops(gu2))]:

gu3:=[seq(MedianTH(i),i=1..N)]:

print(`The first`, N, `terms of the sequence of (truncated) average of the r.v. "total height" defined on rooted labeled trees is`):
print(``):
lprint(gu1):
print(``):

print(`The first`, N, `terms of the sequence of the modes (most likely total height) of the r.v. "total height" defined on rooted labeled trees is`):
print(``):
lprint(gu2a):
print(``):

print(`The first`, N, `terms of the sequence of the (truncated) medians of the r.v. "total height" defined on rooted labeled trees is`):
print(``):
lprint(gu3):
print(``):

print(`Finally, the actual most likely total-heights for i from 1 to`, N, ` are `):
print(``):
lprint(gu2b):
print(``):



end: