######################################################################
##TREES.txt: Save this file as TREES.txt                             #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read TREES.txt<Enter>                              #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Andrew Lohr and Doron Zeilberger, Rutgers University    #
#zeilberg at math dot rutgers dot edu                                #
######################################################################
 
#Created: Dec. 2016-Jan. 2017
 
print(`Created:  Dec. 2016/Jan. 2017`):
print(` This is  TREES.txt `):
print(`It is one of the packages that accompany the article `):
print(` The Limiting Distibutions of the Total Height On Families of Trees`):
print(`by Andrew Lohr and Doron Zeilberger`):

print(`and also available from Lohr's Zeilberger's websites`):
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  procedures from Findrec type ezraF();, for help with`):
 print(`a specific procedure, type ezraF(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

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

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(`---------------------------------------`):

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

Digits:=30:

with(combinat):

ezraBT:=proc()

if args=NULL then
 print(` The Binary Trees procedures are: AsyCn, LimMomsBT, MomBT, MomBTasy, NiBT, NiGFbt, PaperBT  `):
 print(``):

else
ezra(args):
fi:

end:

ezraS:=proc()

if args=NULL then
 print(` The STORY procedures are: Paper1, Paper2  `):
 print(``):

else
ezra(args):
fi:

end:

ezra1:=proc()

if args=NULL then
 print(` The supporting procedures are: AlphaSeq, CheckEqs, CheckJij, Eqs, JijE, Nigzeret, NumSeq, TestNigzeret, Zinn `):
 print(``):

else
ezra(args):
fi:

end:

ezra:=proc()

if args=NULL then
 print(`The main procedures are: Alpha, AlphaF, AlphaFF, AlphaFFw, AlphaFFwNor, Jij, NiGF, NiSeq, NiSeqF, NiSeqFF, NiSeqFFw, PolSeq `):
 print(` `):


elif nops([args])=1 and op(1,[args])=Alpha then
print(`Alpha(P,R,N,k): inputs a polynomial P in the variable R, a variable y, and positive integers N and k, finds`):
print(`the first N terms of the Alpha Seq of the polynomials J_i(y), that are the coefficiens of x^i in the`):
print(`solution of the functional equation`):
print(`J(x,y)=x*P(J(x*y,y)); Done DIRECTLY, very slow! Try`):
print(`Alpha(1+R+R^2,R,20,6);`):

elif nops([args])=1 and op(1,[args])=AlphaF then
print(`AlphaF(P,R,N,k):  A faster version of Alpha(P,R,N,k) (q.v.) using NiSeqF(P,R,N,r) (q.v) for f from 0 to r.`):
print(`Try: `):
print(` AlphaF(1+R+R^2,R,20,4); `):

elif nops([args])=1 and op(1,[args])=AlphaFF then
print(`AlphaFF(P,R,N,k,MaxC):  A yet faster version of AlphaF(P,R,N,k) (q.v.) using recurrences of complexity <=MaxC.`):
print(`Try: `):
print(` AlphaFF(1+R+R^2,R,200,4,16); `):

elif nops([args])=1 and op(1,[args])=AlphaFFw then
print(`AlphaFFw(P,R,N0,N1,k,MaxC):  If N0<N1 are positive integers, the porition of AlphaFF(P,R,N1,k,MaxC):  from n=N0 to n=N1`):
print(`Note that N1-N0 must be at least Max C .`):
print(`Try: `):
print(` AlphaFFw(1+R+R^2,R,2000,2017,4,16); `):


elif nops([args])=1 and op(1,[args])=AlphaFFwNor then
print(`AlphaFFwNor(P,R,N0,N1,k,MaxC):  Similar to AlphaFFw(P,R,N0,N1,k,MaxC) (q.v.), `):
print(`Returns a list of length k`):
print(` whose (1) first entry is the sequence of averages divided by n^(3/2)  for n from N0 to N1`):
print(` (2) second entry is the sequence of coefficient of variations for n from N0 to N1 (s.d. divided by average)` ):
print(`The 3-th through k-th entries are the sequences of alpha-coefficients`):
print(`All these sequenes converge to constants, and the last entries provide appxroximations to the limits`):
print(`Try `):
print(` AlphaFFwNor(1+R+R^2,R,2000,2020,4,16); `):


elif nops([args])=1 and op(1,[args])=AlphaSeq then
print(`AlphaSeq(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`):
prin (` 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(`AlphaSeq((1+x)^100,x,4);`):

elif nops([args])=1 and op(1,[args])=AsyCn then
print(`AsyCn(n,K): the asympototic expansion of binomial(2*n-2,n-1)/(4^n*n^(3/2)) up to order K. Try:`):
print(`AsyCn(n,5);`):

elif nops([args])=1 and op(1,[args])=CheckEqs then
print(`CheckEqs(P,R,m,N): Emprically checks procedure  Eqs(P,R,J,x,m) (q.v.) up to N terms`):
print(` Try: `):
print(`CheckEqs(1+R+R^2,R,2,20);`):

elif nops([args])=1 and op(1,[args])=CheckJij then
print(`CheckJij(P,R,i,j,N):  Checks procedure Jij(P,R,i,j,x,R0) (q.v.) up to N terms.`):
print(` Try: `):
print(` CheckJij(R^2+R+1,R,2,2,30); `):

elif nops([args])=1 and op(1,[args])=Eqs then
print(`Eqs(P,R,J,x,m): the coefficients of t^0 through t^m of J(x,t)-x*P(J(x+x*t,t)) in terms of J[i,r](x), where J[i,r](x) is the -rth derivative`):
print(` of J[r](x),ad J[r](x) is the coeff. of t^r in J(x,t). Try:`):
print(`Eqs(1+R+R^2,R,J,x,2);`):


elif nops([args])=1 and op(1,[args])=JijE then
print(`JijE(P,R,i,j,x,N): inputs a polynomial P in the variable R, and non-neg. integers i and j, and symbols x and R0,`):
print(`outputs the first N terms (i.e. up to x^N) of the formal power series  diff(J_i(x),x$j)`):
print(`is the coeff. of t^i/i! in the solution of the functional equation`):
print(`J(x,t)= x*P(J(x+x*t,t)). Try:`):
print(`JijE(R^2+R+1,R,2,2,x,20);`):

elif nops([args])=1 and op(1,[args])=Nigzeret then
print(`Nigzeret(P,R,x,R0,j): inputs a polynomial P in the variable R, a variable x, and a symbol R0, as well as a positive integer j`):
print(`returns the rational function in R0,x that represents the j-th derivative of the algebaric formal power series`):
print(`defined by R0(x)=x*P(R0(x)). Try:`):
print(`Nigzeret(R^2+R+1,R,x,R0,1);`):


elif nops([args])=1 and op(1,[args])=Jij then
print(`Jij(P,R,i,j,x,R0): inputs a polynomial P in the variable R, and non-neg. integers i and j, and symbols x and R0,`):
print(`outputs a rational function in R0 and x that equals the j-th derivative w.r.t. to x of J_i(x), where J_i(x)`):
print(`is the coeff. of t^i/i! in the solution of the functional equation`):
print(`J(x,t)= x*P(J(x+x*t,t)). Try:`):
print(`Jij(R^2+R+1,R,2,2,x,R0);`):


elif nops([args])=1 and op(1,[args])=LimMomsBT then
print(`LimMomsBT(N): inputs a positive integer N, and outputs a list of numbers, of length N whose`):
print(`first component is the limit of of the average of the total height of complete binary trees with n leaves, divided by n^(3/2),`):
print(`The second component is the limit of the coeff. of variation and the third through the N componets are the`):
print(`limits of the scaled moments, aka as alpha coefficients. Try:`):
print(`LimMomsBT(4);`):

elif nops([args])=1 and op(1,[args])=MomBT then
print(`MomBT(i,n): Inputs a positive integer i, and outputs an explicit expression, in n, for the i-th moment about the mean`):
print(` (unless i=0 when it gives you the actual enumeration, and i=1 when it gives you the average) of the random variable `):
print(` "sum of distances from the root taken over all vertices" defined over all complete binary trees with n leaves (and hence 2n-1 vertices) `):
print(` Try: `):
print(` MomBT(2,n); `):

elif nops([args])=1 and op(1,[args])=MomBTasy then
print(`MomBTasy(i,m,K): Inputs a positive integer i, and outputs an ASYMPTOTIC expression, to order K,`):
print(` in m (that stands for sqrt(n)), for the i-th moment about the mean`):
print(` (unless i=0 when it gives you the actual enumeration, and i=1 when it gives you the average) of the random variable `):
print(` "sum of distances from the root taken over all vertices" defined over all complete binary trees with n leaves (and hence 2n-1 vertices) `):
print(` Try: `):
print(` MomBTasy(1,m,6); `):

elif nops([args])=1 and op(1,[args])=NiBT then
print(`NiBT(i,n): an explicit formula for the sum of the i-th powers of the total heights over all binary trees with 2n-1 vertices (and hence n leaves)`):
print(`w here i is a specfic positive integer, and n is a symbol (denoting a postitive integer). Try: `):
print(`NiBT(2,n); `):


elif nops([args])=1 and op(1,[args])=NiGF then
print(`NiGF(P,R,i,x,R0): The explicit expression in terms of R0 for the generating functon of NiSeq(P,R,N,i) (q.v.)`):
print(`Recall that R0 is the formal power series solution  of the equation R0(x)=x*P(R0(x))`):
print(`Try: `):
print(` NiGF(R^2+1,R,2,x,R0); `):

elif nops([args])=1 and op(1,[args])=NiGFbt then
print(`NiGFbt(i,A): expresses NiGF(1+R^2,R,i,x,R0) in terms of A=sqrt(1-4*x^2). Try:`):
print(`NiGFbt(2,A);`):

elif nops([args])=1 and op(1,[args])=NiSeq then
print(`NiSeq(P,R,N,r) : the first N terms of the numerical sequence (yd/dy)^r P_n(y) where P_n(y) is the n-th term of PolSeq(P,R,y,n) (q.v.)`):
print(`Try:`):
print(`NiSeq(exp(R),R,10,2);`):
print(`NiSeq(1+R^2,R,10,1);`):
print(`NiSeqSeq(1/(1-R),R,10,1);`):

elif nops([args])=1 and op(1,[args])=NiSeqF then
print(`NiSeqF(P,R,N,r) : the first N terms of the numerical sequence (yd/dy)^r P_n(y) where P_n(y) is the n-th term of PolSeq(P,R,y,n) (q.v.)`):
print(`done via the explict generating function given by NiGF(P,R,r,x,R0)`):
print(`Try: `):
print(` NiSeqF(1+R^2,R,20,1); `):
print(` NiSeqF(1+R+R^2,R,20,1); `):

elif nops([args])=1 and op(1,[args])=NiSeqFF then
print(`NiSeqFF(P,R,N0,r,MaxC) : the first N0 terms of the numerical sequence (yd/dy)^r P_n(y) where P_n(y) is the n-th term of PolSeq(P,R,y,n) (q.v.)`):
print(`done really fast, first via the explict generating function given by NiGF(P,R,r,x,R0) and then via the guessed recurrence`):
print(`of complexity MaxC. If no recurrence is found, it returns FAIL, and you should make MaxC larger`):
print(`Try: `):
print(`NiSeqFF(1+R+R^2,R,500,1,10);`):

elif nops([args])=1 and op(1,[args])=NiSeqFFw then
print(`NiSeqFFw(P,R,N0a,N0b, r,MaxC) : the first N0a through N0b terms of the numerical sequence (yd/dy)^r P_n(y)`):
print(` where P_n(y) is the n-th term of PolSeq(P,R,y,n) (q.v.)`):
print(`done really fast, first via the explict generating function given by NiGF(P,R,r,x,R0) and then via the guessed recurrence`):
print(`of complexity MaxC. If no recurrence is found, it returns FAIL, and you should  make MaxC larger`):
print(`Note that N0b-N0a must be at least Max C. `):
print(`Try: `):
print(`NiSeqFFw(1+R+R^2,R,500,511,1,10);`):

elif nops([args])=1 and op(1,[args])=NumSeq then
print(`NumSeq(P,R,N): inputs a polynomial (or even a function) P of R, outputs the list of numbers`):
print(`of length N, whose i-th entry is the coeff. of x^i in the solution of the algebraic equation`):
print(`R(x)=x*P(R(x)); Try: `):
print(` NumSeq(1+R^2,R,10); `):
print(` NumSeq(1/(1-R),R,10); `):


elif nops([args])=1 and op(1,[args])=Paper1 then
print(`Paper1(DegSet,k,MaxC,N0,N1,N2): input a set, DegSet, of positive integers, and symbols R,x,y,R0, `):
print(` and posivie integers k and MaxC, and large positive integers N0 and N1,  and outputs an article about `):
print(` the statistical behavior of the random variable "sum of heights of all vertices" defined in the sample space of `):
print(` of all n-vertex ordered rooted trees where each vertex is either a leaf (no children) or must have a number of children `):
print(` in DegSet. It gives you information about all the moments up to the k-th. MaxC is the maximum complexity of the `):
print(` recurrence looked for the numerators of the straight moments. `):
print(` For example, for rooted ordered trees where each vertex must have 0,1,or 2, children, type `):
print(` Paper1({1,2},4,20,20,2000,2020); `):


elif nops([args])=1 and op(1,[args])=Paper2 then
print(`Paper2(DegSet,k,N0): Like Paper1(DegSet,k,MaxC,N0,N1,N2) but only generating functions.`):
print(` For example, for rooted ordered trees where each vertex must have 0,1,or 2, children, type `):
print(` Paper2({1,2},4,20); `):

elif nops([args])=1 and op(1,[args])=PaperBT then
print(`PaperBT(N): input a positive integer N and outputs an article about the explicit expressions, as well as asymptotic`):
print(`for the average, variance, and the third through the N-th moment, as well as the of the scaled limits`):
print(`for the random variable: total height defined over complete binary trees with n leaves. Try:`):
print(`PaperBT(4):`):

elif nops([args])=1 and op(1,[args])=PolSeq then
print(`PolSeq(P,R,y,N): inputs a polynomial (or even a function) P of R, outputs the list of polynomials in y`):
print(`of length N, whose i-th entry is the coeff. of x^i in the solution of the functional equation`):
print(`J(x,y)=x*P(J(x*y,y)); Try:`):
print(`PolSeq(exp(R),R,y,10);`):
print(`PolSeq(1+R^2,R,y,10);`):
print(`PolSeq(1/(1-R),R,y,10);`):
print(``):


elif nops([args])=1 and op(1,[args])=TestNigzeret then
print(`TestNigzeret(P,R,J,N): tests procedure Nigzeret(P,R,x,R0,j,N) (q.v.) for  j from 1 to J up to N terms`):
print(`Try: `):
print(` TestNigzeret(R^2+R+1,R,4,40); `):

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:


####From FindRec
######################################################################
## FindRec: Save this file as FindRec To use it, stay in the         #
## same directory, get into Maple (by typing: maple <Enter> )        #
## and then type:  read Findrec <Enter>                              #
## Then follow the instructions given there                          #
##                                                                   #
## Written by Doron Zeilberger, Rutgers University ,                 #
##  zeilberg@math.rutgers.edu.                                       # 
######################################################################

 
ezraF:=proc()
if args=NULL then

print(` FindRec: A Maple package for empirically guessing partial recurrence`):
print(`equations satisfied by Discrete Functions of ONE Variable`):
print():
print(`For help with a specific procedure, type "ezra(procedure_name);"`):
print(`Contains procedures:  `):
print(`  findrec, Findrec, FindrecF, SeqFromRec `):
print():

elif nargs=1 and args[1]=findrec then
print(`findrec(f,DEGREE,ORDER,n,N): guesses a recurrence operator annihilating`):
print(`the sequence f of degree DEGREE and order ORDER.`):
print(`For example, try: findrec([seq(i,i=1..10)],0,2,n,N);`):

elif nargs=1 and args[1]=Findrec then
print(`Findrec(f,n,N,MaxC): Given a list f tries to find a linear recurrence equation with`):
print(`poly coffs. ope(n,N), where n is the discrete variable, and N is the shift operator `):
print(`of maximum DEGREE+ORDER<=MaxC`):
print(`e.g. try Findrec([1,1,2,3,5,8,13,21,34,55,89],n,N,2);`):

elif nargs=1 and args[1]=FindrecF then
print(`FindrecF(f,n,N): Given a function f of a single variable tries to find a linear recurrence equation with`):
print(`poly coffs. .g. try FindrecF(i->i!,n,N);`):


elif nargs=1 and args[1]=SeqFromRec then
print(`SeqFromRec(ope,n,N,Ini,K): Given the first L-1`):
print(`terms of the sequence Ini=[f(1), ..., f(L-1)]`):
print(`satisfied by the recurrence ope(n,N)f(n)=0`):
print(`extends it to the first K values`):
print(`For example, try:`):
print(`SeqFromRec(N-n-1,n,N,[1],10);`):


fi:
 
end:


###Findrec 
#findrec(f,DEGREE,ORDER,n,N): guesses a recurrence operator annihilating
#the sequence f of degree DEGREE and order ORDER
#For example, try: findrec([seq(i,i=1..10)],0,2,n,N);
findrecVerbose:=proc(f,DEGREE,ORDER,n,N)
local ope,var,eq,i,j,n0,kv,var1,eq1,mu,a:
if (1+DEGREE)*(1+ORDER)+3+ORDER>nops(f) then
ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE):
fi:
ope:=0:
var:={}:
 
for i from 0 to ORDER do
 for j from 0 to DEGREE do
  ope:=ope+a[i,j]*n^j*N^i:
  var:=var union {a[i,j]}:
 od:
od:
    
eq:={}:
 
for n0 from 1 to (1+DEGREE)*(1+ORDER)+2 do
  eq1:=0:
 
  for i from 0 to ORDER do
     eq1:=eq1+subs(n=n0,coeff(ope,N,i))*op(n0+i,f):
  od:
 
   eq:= eq union {eq1}:
od:
 
var1:=solve(eq,var):
 
kv:={}:
 
for i from 1 to nops(var1) do
 mu:=op(i,var1):
 
if op(1,mu)=op(2,mu) then
   kv:= kv union {op(1,mu)}:
 fi:
 
od:

ope:=subs(var1,ope):

if ope=0 then
  RETURN(FAIL):
fi:

ope:={seq(coeff(expand(ope),kv[i],1),i=1..nops(kv))} minus {0}:

if nops(ope)>1 then
print(`There is some slack, there are `, nops(ope)):
print(ope):
RETURN(Yafe(ope[1],N)[2]):
elif nops(ope)=1 then
RETURN(Yafe(ope[1],N)[2]):
else
 RETURN(FAIL):
fi:

end:
 
#findrec(f,DEGREE,ORDER,n,N): guesses a recurrence operator annihilating
#the sequence f of degree DEGREE and order ORDER
#For example, try: findrec([seq(i,i=1..10)],0,2,n,N);
findrec:=proc(f,DEGREE,ORDER,n,N)
local ope,var,eq,i,j,n0,kv,var1,eq1,mu,a:
option remember:


if not findrecEx(f,DEGREE,ORDER,ithprime(20)) then
 RETURN(FAIL):
fi:

if not findrecEx(f,DEGREE,ORDER,ithprime(40)) then
 RETURN(FAIL):
fi:

if not findrecEx(f,DEGREE,ORDER,ithprime(80)) then
 RETURN(FAIL):
fi:


if (1+DEGREE)*(1+ORDER)+5+ORDER>nops(f) then
ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE):
fi:
ope:=0:
var:={}:
 
for i from 0 to ORDER do
 for j from 0 to DEGREE do
  ope:=ope+a[i,j]*n^j*N^i:
  var:=var union {a[i,j]}:
 od:
od:
    
eq:={}:
 
for n0 from 1 to (1+DEGREE)*(1+ORDER)+4 do
  eq1:=0:
 
  for i from 0 to ORDER do
     eq1:=eq1+subs(n=n0,coeff(ope,N,i))*op(n0+i,f):
  od:
 
   eq:= eq union {eq1}:
od:
 
var1:=solve(eq,var):
 
kv:={}:
 
for i from 1 to nops(var1) do
 mu:=op(i,var1):
 
if op(1,mu)=op(2,mu) then
   kv:= kv union {op(1,mu)}:
 fi:
 
od:

ope:=subs(var1,ope):

if ope=0 then
  RETURN(FAIL):
fi:

ope:={seq(coeff(expand(ope),kv[i],1),i=1..nops(kv))} minus {0}:

if nops(ope)>1 then
RETURN(Yafe(ope[1],N)[2]):
elif nops(ope)=1 then
RETURN(Yafe(ope[1],N)[2]):
else
 RETURN(FAIL):
fi:

end:
 


Yafe:=proc(ope,N) local i,ope1,coe1,L: 
if ope=0 then
 RETURN(1,0):
fi:
ope1:=expand(ope):
L:=degree(ope1,N):
coe1:=coeff(ope1,N,L):
ope1:=normal(ope1/coe1):
ope1:=normal(ope1):
ope1:=
convert(
[seq(factor(coeff(ope1,N,i))*N^i,i=ldegree(ope1,N)..degree(ope1,N))],`+`):
factor(coe1),ope1:
end:


#Findrec(f,n,N,MaxC): Given a list f tries to find a linear recurrence equation with
#poly coffs.
#of maximum DEGREE+ORDER<=MaxC
#e.g. try Findrec([1,1,2,3,5,8,13,21,34,55,89],n,N,2);
Findrec:=proc(f,n,N,MaxC)
local DEGREE, ORDER,ope,L:

for L from 0 to MaxC do
for ORDER from 1 to L do
 DEGREE:=L-ORDER:
if (2+DEGREE)*(1+ORDER)+4>=nops(f) then
 print(`Insufficient data for degree`, DEGREE, `and order `,ORDER):
 RETURN(FAIL):
fi:

 ope:=findrec([op(1..(2+DEGREE)*(1+ORDER)+4,f)],DEGREE,ORDER,n,N):
     if ope<>FAIL then
       RETURN(ope):
      fi:
 od:
od:
FAIL:

end:





 
#SeqFromRec(ope,n,N,Ini,K): Given the first L-1
#terms of the sequence Ini=[f(1), ..., f(L-1)]
#satisfied by the recurrence ope(n,N)f(n)=0
#extends it to the first K values
SeqFromRec:=proc(ope,n,N,Ini,K)
local ope1,gu,L,n1,j1:
ope1:=Yafe(ope,N)[2]:
L:=degree(ope1,N):
if nops(Ini)<>L then
 ERROR(`Ini should be of length`, L):
fi:

ope1:=expand(subs(n=n-L,ope1)/N^L):

gu:=Ini:

for n1 from nops(Ini)+1 to K do
gu:=[op(gu), -add(gu[nops(gu)+1-j1]*subs(n=n1,coeff(ope1,N,-j1)),
j1=1..L)]:
od:

gu:

end:

#SeqFromRecW(ope,n,N,Ini,K1,K2): Given the first L-1
#terms of the sequence Ini=[f(1), ..., f(L-1)]
#satisfied by the recurrence ope(n,N)f(n)=0
#finds the values from K1 to K2
SeqFromRecW:=proc(ope,n,N,Ini,K1,K2)
local ope1,gu,L,n1,j1,w:
w:=K2-K1+1:
ope1:=Yafe(ope,N)[2]:
L:=degree(ope1,N):
if nops(Ini)<>L then
 ERROR(`Ini should be of length`, L):
fi:

ope1:=expand(subs(n=n-L,ope1)/N^L):

gu:=Ini:

for n1 from nops(Ini)+1 to nops(Ini)+w do
gu:=[op(gu), -add(gu[nops(gu)+1-j1]*subs(n=n1,coeff(ope1,N,-j1)),
j1=1..L)]:
od:

gu:=[op(nops(gu)-w+1..nops(gu),gu)]:


for n1 from  nops(Ini)+w+1 to K2  do
gu:=[op(gu), -add(gu[nops(gu)+1-j1]*subs(n=n1,coeff(ope1,N,-j1)),
j1=1..L)]:
gu:=[op(2..nops(gu),gu)]:
od:
gu:

end:
 
#end Findrec


with(linalg):

#findrecEx(f,DEGREE,ORDER,m1): Explores whether thre
#is a good chance that there is a recurrence of degree DEGREE
#and order ORDER, using the prime m1
#For example, try: findrecEx([seq(i,i=1..10)],0,2,n,N,1003);
findrecEx:=proc(f,DEGREE,ORDER,m1)
local ope,var,eq,i,j,n0,eq1,a,A1,
D1,E1,Eq,Var,f1,n,N:
option remember:
f1:=f mod m1:
if (1+DEGREE)*(1+ORDER)+5+ORDER>nops(f) then
ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE):
fi:
ope:=0:
var:={}:
 
for i from 0 to ORDER do
 for j from 0 to DEGREE do
  ope:=ope+a[i,j]*n^j*N^i:
  var:=var union {a[i,j]}:
 od:
od:
    
eq:={}:
 
for n0 from 1 to (1+DEGREE)*(1+ORDER)+4 do
  eq1:=0:
 
  for i from 0 to ORDER do
     eq1:=eq1+subs(n=n0,coeff(ope,N,i))*op(n0+i,f1) mod m1:
  od:
 
   eq:= eq union {eq1}:
od:


Eq:= convert(eq,list):
Var:= convert(var,list):


D1:=nops(Var):
E1:=nops(Eq):
if E1<D1 then
  RETURN(true):
fi:

A1:=matrix(D1,D1):

for i from 1 to D1-1 do
for j from 1 to D1 do
  A1[i,j]:=coeff(Eq[i],Var[j]):
od:
od:

for j from 1 to nops(Var) do
  A1[D1,j]:=coeff(Eq[D1],Var[j]):
od:

if det(A1) mod m1 <>0  then
 RETURN(false):
fi:

if E1-D1>=1 then
 for j from 1 to nops(Var) do
  A1[D1,j]:=coeff(Eq[D1+1],Var[j]):
 od:

if det(A1) mod m1 <>0  then
 RETURN(false):
fi:
fi:

if E1-D1>=2 then
 for j from 1 to nops(Var) do
  A1[D1,j]:=coeff(Eq[D1+2],Var[j]):
 od:

if det(A1) mod m1 <>0  then
 RETURN(false):
fi:
fi:

true:

end:
####End From FindRec

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,50 ):
mu:=evalf(sqrt(op(n1+1,resh)/op(n1-1,resh))*exp(-(s1+s2)/((s1-s2)*s1)),50):
evalf([theta,mu],10):
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:

#AlphaSeq(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:
#AlphaSeq(((1+x)/2)^100,x,4);
AlphaSeq:=proc(f,x,N) local gu,i:
 gu:=AveAndMoms(f,x,N):
 if gu=FAIL then
  RETURN(gu):
 fi:

if gu[2]=0 then
  RETURN([gu[1]]):
fi:

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

end:


#PolSeq(P,R,y,N): inputs a polynomial (or even a function) P of R, outputs the list of polynomials in y
#of length N, whose i-th entry is the coeff. of x^i in the solution of the functional equation
#J(x,y)=x*P(J(x*y,y)); Try
#PolSeq(exp(R),R,y,10);
#PolSeq(1+R^2,R,y,10);
#PolSeq(1/(1-R),R,y,10);
PolSeq:=proc(P,R,y,N) local x,gu,gu1,gu2,i:
option remember:
gu:=x:

gu1:=x*subs(R=subs(x=x*y,gu),P):
gu1:=taylor(gu1,x=0,N+1):
gu1:=add(expand(coeff(gu1,x,i))*x^i,i=1..N):

while gu1<>gu do
gu2:=x*subs(R=subs(x=x*y,gu1),P):
gu2:=taylor(gu2,x=0,N+1):
gu2:=add(expand(coeff(gu2,x,i))*x^i,i=1..N):
gu:=gu1:
gu1:=gu2:
od:
[seq(coeff(gu,x,i),i=1..N)]:
end:

#NumSeq(P,R,N): inputs a polynomial (or even a function) P of R, outputs the list of numbers
#of length N, whose i-th entry is the coeff. of x^i in the solution of the algebraic equation
#R(x)=x*P(R(x)); Try
#NumSeq(exp(R),R,10);
#NumSeq(1+R^2,R,10);
#NumSeq(1/(1-R),R,10);
NumSeq:=proc(P,R,N) local x,gu,gu1,gu2,i:
option remember:
gu:=x:

gu1:=x*subs(R=gu,P):
gu1:=taylor(gu1,x=0,N+1):
gu1:=add(expand(coeff(gu1,x,i))*x^i,i=1..N):

while gu1<>gu do
gu2:=x*subs(R=gu1,P):
gu2:=taylor(gu2,x=0,N+1):
gu2:=add(expand(coeff(gu2,x,i))*x^i,i=1..N):
gu:=gu1:
gu1:=gu2:
od:
[seq(coeff(gu,x,i),i=1..N)]:
end:








#NiSeq(P,R,N,r) : the first N terms of the numerical sequence (yd/dy)^r P_n(y) where P_n(y) is the n-th term of PolSeq(P,R,y,n) (q.v.)
#Try:
#NiSeq(exp(R),R,y,10,2);
#NiSeq(1+R^2,R,1);
#NiSeq(1/(1-R),R,1);
NiSeq:=proc(P,R,N,r)  local y,gu,i,i1:
gu:=PolSeq(P,R,y,N):
for i from 1 to r do
 gu:=expand([seq(y*diff(gu[i1],y),i1=1..nops(gu))]):
od:
subs(y=1,gu):
end:




#Nigzeret(P,R,x,R0,j): inputs a polynomial P in the variable R, a variable x, and a symbol R0, as well as a positive integer j
#returns the rational function in R0,x that represents the j-th derivative of the algebaric formal power series
#defined by R0(x)=x*P(R0(x)). Try:
#Nigzeret(R^2+R+1,R,x,R0,1);
Nigzeret:=proc(P,R,x,R0,j) local A,gu,R1,mu:
option remember:

if j=0 then
 RETURN(R0):
fi:

if j=1 then

gu:=A(x)-x*subs(R=A(x),P):
gu:=diff(gu,x):
gu:=subs({diff(A(x),x)=R1, A(x)=R0},gu):
R1:=solve(gu,R1):
R1:=rem(numer(R1),R0-x*subs(R=R0,P),R0)/rem(denom(R1),R0-x*subs(R=R0,P),R0):
return(R1):
fi:

mu:=Nigzeret(P,R,x,R0,1):

gu:=Nigzeret(P,R,x,R0,j-1):
gu:=subs(R0=A(x),gu):
gu:=diff(gu,x):
gu:=normal(subs({diff(A(x),x)=mu,A(x)=R0}, gu)):
gu:=normal(rem(numer(gu),R0-x*subs(R=R0,P),R0)/rem(denom(gu),R0-x*subs(R=R0,P),R0)):

RETURN(gu):

end:

#TestNigzeret(P,R,J,N): tests procedure Nigzeret(P,R,x,R0,j,N) (q.v.) for  j from 1 to J up to N terms
#Try
#TestNigzeret(R^2+R+1,R,4,40);

TestNigzeret:=proc(P,R,J,N) local gu0,x,R0,i,j,mu:

 gu0:=PolSeq(P,R,1,N+J+2):
 gu0:=add(gu0[i]*x^i,i=1..N+J+2):


  for j from 1 to J do
   mu:=Nigzeret(P,R,x,R0,j):
   mu:=subs(R0=gu0,mu):
   mu:=expand(diff(gu0,x$j)-mu):
   mu:=taylor(mu,x,N+J):
  if {seq(expand(coeff(mu,x,i)),i=1..N)}<>{0} then
     RETURN(false):
  fi:
 od:

true:

end:


#Eqs(P,R,J,x,m): the coefficients of t^0 through t^m of J(x,t)-x*P(J(x+x*t,t)) in terms of J[i,r](x), where J[i,r](x) is the -rth derivative
# of J[r](x),ad J[r](x) is the coeff. of t^r in J(x,t). Try:
#Eqs(1+R+R^2,R,J,x,2);
Eqs:=proc(P,R,J,x,m) local gu,i,t,r,gu1:
gu:=add(J[i,0]*t^i/i!,i=0..m):
gu1:=add( add(J[i,r]/r!*(x*t)^r,r=0..m)*t^i/i!,i=0..m):
gu:=expand(gu-x*subs(R=gu1,P)):
[seq(expand(coeff(gu,t,i)),i=0..m)]:
end:


#CheckEqs(P,R,m,N): Emprically checks procedure  Eqs(P,R,J,x,m) (q.v.) up to N terms
#Try:
#CheckEqs(1+R+R^2,R,2,20);
CheckEqs:=proc(P,R,m,N) local gu,J,x,i,j:
gu:=Eqs(P,R,J,x,m):

for i from 0 to m do
 for j from 0 to m do
  gu:=expand(subs(J[i,j]=JijE(P,R,i,j,x,N),gu)):
 od:
od:

evalb({seq(seq(coeff(gu[i1],x,j),j=0..N),i1=1..nops(gu))}={0}):
 
end:


#Jij(P,R,i,j,x,R0): inputs a polynomial P in the variable R, and non-neg. integers i and j, and symbols x and R0,
#outputs a rational function in R0 and x that equals the j-th derivative w.r.t. to x of J_i(x), where J_i(x)
#is the coeff. of t^i/i! in the solution of the functional equation
#J(x,t)= x*P(J(x+x*t,t)). Try:
#Jij(R^2+R+1,R,2,2,x,R0);
Jij:=proc(P,R,i,j,x,R0) local gu,mu,A,i1,J,j1:
option remember:
 if i=0 then
  RETURN(Nigzeret(P,R,x,R0,j)):
 fi:

mu:=Nigzeret(P,R,x,R0,1):

if j>0 then
 gu:=normal(Jij(P,R,i,j-1,x,R0)):
 gu:=subs(R0=A(x),gu):
 gu:=normal(diff(gu,x)):
 gu:=subs({A(x)=R0,diff(A(x),x)=mu},gu):
 gu:=normal(gu):
gu:=normal(rem(numer(gu),R0-x*subs(R=R0,P),R0)/rem(denom(gu),R0-x*subs(R=R0,P),R0)):
 RETURN(gu):
fi:

gu:=Eqs(P,R,J,x,i)[i+1]:
gu:=normal(solve(gu,J[i,0])):

for i1 from 0 to i-1 do
for j1 from 0 to i+2 do
 gu:=subs(J[i1,j1]=Jij(P,R,i1,j1,x,R0),gu):
od:
od:
gu:=normal(gu):
gu:=normal(rem(numer(gu),R0-x*subs(R=R0,P),R0)/rem(denom(gu),R0-x*subs(R=R0,P),R0)):
 RETURN(gu):

end:
  

 

#CheckJij(P,R,i,j,N):  Checks procedure Jij(P,R,i,j,x,R0) (q.v.) up to N terms.
#Try:
#CheckJij(R^2+R+1,R,2,2,30);
CheckJij:=proc(P,R,i,j,N) local gu,x,R0,mu,y,mu0,i1,i2:

gu:=Jij(P,R,i,j,x,R0):
mu:=PolSeq(P,R,y,N+3*j+3*i+10):
mu0:=subs(y=1,mu):
mu0:=add(mu0[i2]*x^i2,i2=1..N+3*j+3*i+10):

for i1 from 1 to i do
mu:=[seq(diff(mu[i2],y),i2=1..nops(mu))]:
od:

mu:=subs(y=1,mu):


mu:=add(mu[i1]*x^i1,i1=1..N+3*j+3*i+10):

if j>0 then
 mu:=diff(mu,x$j):
fi:

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

gu:=normal(subs(R0=mu0,gu)):

gu:=taylor(gu,x=0,N+1):
gu:=[seq(coeff(gu,x,i2),i2=1..N)]:

if gu=mu then
 print(gu,mu):
 RETURN(true):
else
 print(gu,mu):
 RETURN(false):
fi:

end:





#JijE(P,R,i,j,x,N): inputs a polynomial P in the variable R, and non-neg. integers i and j, and symbols x and R0,
#outputs the first N terms (i.e. up to x^N) of the formal power series  diff(J_i(x),x$j)
#is the coeff. of t^i/i! in the solution of the functional equation
#J(x,t)= x*P(J(x+x*t,t)). Try:
#Jij(R^2+R+1,R,2,2,x,R0);
JijE:=proc(P,R,i,j,x,N) local gu,i1,i2,y:
gu:=PolSeq(P,R,y,N+j):

for i1 from 1 to i do
 gu:=[seq(diff(gu[i2],y),i2=1..N+j)]:
od:

gu:=subs(y=1,gu):
gu:=add(gu[i1]*x^i1,i1=1..N+j):

if j>0 then
 gu:=diff(gu,x$j):
fi:

gu:
end:

 





#NiGF(P,R,i,x,R0): The explicit expression in terms of R0 for the generating functon of NiSeq(P,R,N,i) (q.v.)
#Recall that R0 is the formal power series solution  of the equation R0(x)=x*P(R0(x))
#Try:
#NiGF(R^2+1,R,2,x,R0);
NiGF:=proc(P,R,i,x,R0) local i1,gu:
option remember:
gu:=normal(add(Jij(P,R,i1,0,x,R0)*stirling2(i,i1),i1=0..i)):
gu:=normal(rem(numer(gu),R0-x*subs(R=R0,P),R0)/rem(denom(gu),R0-x*subs(R=R0,P),R0)):
gu:
end:



#NiSeqF(P,R,N,r) : the first N terms of the numerical sequence (yd/dy)^r P_n(y) where P_n(y) is the n-th term of PolSeq(P,R,y,n) (q.v.)
#done via the explict generating function given by NiGF(P,R,r,x,R0)
#Try:
#NiSeqF(exp(R),R,y,10,2);
#NiSeqF(1+R^2,R,1);
#NiSeqF(1/(1-R),R,1);
NiSeqF:=proc(P,R,N,r)  local mu0,x,R0,i1,gu:
mu0:=NumSeq(P,R,N):
mu0:=add(mu0[i1]*x^i1,i1=1..N):
gu:=NiGF(P,R,r,x,R0):
gu:=subs(R0=mu0,gu):
gu:=taylor(gu,x=0,N+1):
[seq(coeff(gu,x,i1),i1=1..N)]:
end:





#Alpha(P,R,N,k): inputs a polynomial P in the variable R, a variable y, and positive integers N and k, finds
#the first N terms of the Alpha Seq of the polynomials J_i(y), that are the coefficiens of x^i in the
#solution of the functional equation
#J(x,y)=x*P(J(x*y,y)); Try
#Alpha(1+R+R^2,R,20,6);
Alpha:=proc(P,R,N,k) local y,gu,i,ku:
 gu:=PolSeq(P,R,y,N):

if {seq(gu[2*i],i=1..nops(gu)/2)}={0} then
 gu:=[seq(gu[2*i-1],i=1..nops(gu)/2)]:
fi:

if {seq(gu[2*i-1],i=1..nops(gu)/2 )}={0} then
 gu:=[seq(gu[2*i],i=1..nops(gu)/2 )]:
fi:


[seq(AlphaSeq(gu[i],y,k),i=1..nops(gu))]:
end:

#AlphaF(P,R,N,k):  A faster version of Alpha(P,R,N,k) (q.v.) using NiSeqF(P,R,N,r) (q.v) for f from 0 to r.
#Try
#AlphaF(1+R+R^2,R,20,4);
AlphaF:=proc(P,R,N,k) local mu,lu,r,i1,i,ku:

mu[0]:=NiSeqF(P,R,N,0):

for i from 1 to k do
 mu[i]:=NiSeqF(P,R,N,i):
od:


ku:=[]:

for i1 from 1 to N do
 if mu[0][i1]=0 then
  ku:=[op(ku),FAIL]:
 else
  ku:=[op(ku),mu[1][i1]/mu[0][i1]]:
 fi:
od:

lu[1]:=ku:

for i from 2 to k do

ku:=[]:

for i1 from 1 to N do
 if mu[0][i1]=0 then
  ku:=[op(ku),FAIL]:
 else
  ku:=[op(ku),add((-1)^r*binomial(i,r)*mu[1][i1]^r*mu[0][i1]^(i-r-1)*mu[i-r][i1],r=0..i)]:
 fi:

lu[i]:=ku:
od:



 ku:=[]:

for i1 from 1 to N do
 if mu[0][i1]=0 then
  ku:=[op(ku),FAIL]:
 else
  ku:=[op(ku),lu[i][i1]/mu[0][i1]^i ]:
 fi:
od:

 lu[i]:=ku:

od:

lu:=[seq([seq(lu[i][i1],i=1..k)],i1=1..N)]:

ku:=[]:

for i1 from 1 to N do

 if lu[i1][2]=0 or lu[i1][1]=0 then
   ku:=[op(ku),[lu[i1][1]] ]:
else
 ku:=[op(ku), [lu[i1][1],lu[i1][2],seq(lu[i1][i]/lu[i1][2]^(i/2),i=3..k)] ]:
fi:
od:

ku:
end:








#NiSeqFF(P,R,N0,r,MaxC) : the first N0 terms of the numerical sequence (yd/dy)^r P_n(y) where P_n(y) is the n-th term of PolSeq(P,R,y,n) (q.v.)
#done really fast, forst via the explict generating function given by NiGF(P,R,r,x,R0) and then via the guessed recurrence
#of complexity MaxC. If no recurrence is found, it returns FAIL, and you should try to make MaxC larger/
#Try:
#NiSeqFF(1+R+R^2,R,500,1,12);
NiSeqFF:=proc(P,R,N0,r,MaxC)  local gu,L,i,N1,n,N,ope:

N1:=max(seq(seq((2+L-i)*(1+i)+4,i=0..L),L=0..MaxC))+10:

gu:=NiSeqF(P,R,N1,r):

ope:=Findrec(gu,n,N,MaxC):

if ope=FAIL then
  print(` Make `, MaxC, ` bigger `):
  RETURN(FAIL):
fi:

SeqFromRec(ope,n,N,[op(1..degree(ope,N),gu)],N0):

end:

#NiSeqFFw(P,R,N0a,N0b,r,MaxC) : the terms N0a to N0b  of the numerical sequence (yd/dy)^r P_n(y) where P_n(y) is the n-th term of PolSeq(P,R,y,n) (q.v.)
#done really fast, forst via the explict generating function given by NiGF(P,R,r,x,R0) and then via the guessed recurrence
#of complexity MaxC. If no recurrence is found, it returns FAIL, and you should try to make MaxC larger/
#Try:
#NiSeqFFw(1+R+R^2,R,500,510,1,12);
NiSeqFFw:=proc(P,R,N0a,N0b,r,MaxC)  local gu,L,i,N1,n,N,ope:

if N0b-N0a<MaxC then
 print(`N0b-N0a must be at least `, MaxC):
 RETURN(FAIL):
fi:


N1:=max(seq(seq((2+L-i)*(1+i)+4,i=0..L),L=0..MaxC))+10:

gu:=NiSeqF(P,R,N1,r):

ope:=Findrec(gu,n,N,MaxC):

if ope=FAIL then
  print(` Make `, MaxC, ` bigger `):
  RETURN(FAIL):
fi:

SeqFromRecW(ope,n,N,[op(1..degree(ope,N),gu)],N0a,N0b):

end:



#AlphaFF(P,R,N,k,MaxC):  A faster version of AlphaF(P,R,N,k) using recurrences with complexity <=MaxC
#Try
#AlphaFF(1+R+R^2,R,20,4);
AlphaFF:=proc(P,R,N,k,MaxC) local mu,lu,r,i1,i,ku:

mu[0]:=NiSeqFF(P,R,N,0,MaxC):

if mu[0]=FAIL then
  RETURN(FAIL):
fi:

for i from 1 to k do
 mu[i]:=NiSeqFF(P,R,N,i,MaxC):

if mu[i]=FAIL then
  RETURN(FAIL):
fi:

od:

lu[1]:=[seq(mu[1][i]/mu[0][i],i=1..N)]:

for i from 2 to k do
 lu[i]:=[seq(add((-1)^r*binomial(i,r)*mu[1][i1]^r*mu[0][i1]^(i-r-1)*mu[i-r][i1],r=0..i),i1=1..N)]:

 lu[i]:=[seq(lu[i][i1]/mu[0][i1]^i,i1=1..N)]:

od:

lu:=[seq([seq(lu[i][i1],i=1..k)],i1=1..N)]:

ku:=[]:

for i1 from 1 to N do

 if lu[i1][2]=0 or lu[i1][1]=0 then
   ku:=[op(ku),[lu[i1][1]] ]:
else
 ku:=[op(ku), [lu[i1][1],lu[i1][2],seq(lu[i1][i]/lu[i1][2]^(i/2),i=3..k)] ]:
fi:
od:

ku:
end:

#AlphaFFw(P,R,N0,N1,k,MaxC):  When N0<N1 The last portion of AlphaFF(P,R,N1,k,MaxC) starting at n=N0
#Try
#AlphaFFw(1+R+R^2,R,2000,2010,4,16);
AlphaFFw:=proc(P,R,N0,N1,k,MaxC) local mu,lu,r,i1,i,ku,w:

if N1-N0<MaxC then
 print(`N1-N0 must be at least `,MaxC):
 RETURN(FAIL):
fi:

mu[0]:=NiSeqFFw(P,R,N0,N1,0,MaxC):

w:=nops(mu[0]):

if mu[0]=FAIL then
  RETURN(FAIL):
fi:

for i from 1 to k do
 mu[i]:=NiSeqFFw(P,R,N0,N1,i,MaxC):

if mu[i]=FAIL then
  RETURN(FAIL):
fi:

od:

lu[1]:=[seq(mu[1][i]/mu[0][i],i=1..w)]:

for i from 2 to k do
 lu[i]:=[seq(add((-1)^r*binomial(i,r)*mu[1][i1]^r*mu[0][i1]^(i-r-1)*mu[i-r][i1],r=0..i),i1=1..w)]:

 lu[i]:=[seq(lu[i][i1]/mu[0][i1]^i,i1=1..w)]:

od:

lu:=[seq([seq(lu[i][i1],i=1..k)],i1=1..w)]:

ku:=[]:

for i1 from 1 to w do
 ku:=[op(ku), [lu[i1][1],lu[i1][2],seq(lu[i1][i]/lu[i1][2]^(i/2),i=3..k)] ]:

od:


ku:
end:


#AlphaFFwNor(P,R,N0,N1,k,MaxC):  Similar to AlphaFFw(P,R,N0,N1,k,MaxC) (q.v.), 
#Returns a list of length k
#whose (1)first entry is the sequence of averages divided by n^(3/2) 
#(2) second entry is the sequence of coefficient of variations for n from N0 to N1 (s.d. divided by averagge)
#The 3-th through k-th entries are the sequences of alpha-coefficients
#All these sequenes converge to constants, and the last entry is an appxroximation
#Try
#AlphaFFwNor(1+R+R^2,R,2000,2010,4,16);
AlphaFFwNor:=proc(P,R,N0,N1,k,MaxC) local gu,T,gu1,gu1a,gu2,i,i1:

gu:=AlphaFFw(P,R,N0,N1,k,MaxC):

for i from 1 to k do
 T[i]:=evalf(evalf([seq(gu[i1][i],i1=1..nops(gu))],30),20):
od:

gu1:=T[1]:

gu1a:=evalf(evalf([seq(gu1[i1]/(N0+i1-1)^(3/2),i1=1..nops(gu1))],30),20):

gu2:=T[2]:

gu2:=evalf(evalf([seq(sqrt(gu2[i1])/gu1[i1],i1=1..nops(gu2))],30),10):

evalf([gu1a,gu2,seq(T[i],i=3..k)],20):

end:

#Paper1(DegSet,k,MaxC,N0, N1,N2): input a set, DegSet, of positive integers, and symbols R,x,y,R0, 
#and posivie integers k and MaxC, and  a small positive integer N0, and large positive integers N1 and N2,  and outputs an article about
#the statistical behavior of the random variable "sum of heights of all vertices" defined in the sample space of
#of all n-vertex ordered rooted trees where each vertex is either a leaf (no children) or must have a number of children
#in DegSet. It gives you information about all the moments up to the k-th. MaxC is the maximum complexity of the
#recurrence looked for the numerators of the straight moments.
#For example, for rooted ordered trees where each vertex must have 0,1,or 2, children, type
#Paper1({1,2},4,20,20,2000,2020);
Paper1:=proc(DegSet,k,MaxC,N0,N1,N2) local R,x,R0,gu,y,P,i,n,N,R1,t0,i1:

t0:=time():

print(``):
print(`On the Statistical Distribution of the Total Height  on Ordered Trees where each vertex must have a Number of Children in the set`, DegSet union {0}):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`Let T(n) be the set of ordered (rooted) trees on n vertices where each vertex may be a leaf or else must have a number of children in the set`, DegSet):
print(``):
P:=1+add(R0(x)^i, i in DegSet):

print(`Let R0(x) be the ordinary generating function of the enumerating sequence. By generatingfunctionlogy, R0(x) is the unique formal power series`):
print(`satisfying `):
print(``):
print(R0(x)= x*P):
print(``):
print(`The first`, N0, `terms of the enumerating sequence starting at n=1, are `):

gu:=NumSeq(1+add(R1^i,i in DegSet),R1,N0):
print(``):
print(gu):
print(``):
print(`Let H[n] be the random variable "sum of heights (distances from the root) of all vertices", and let P[n](y) be the`):
print(`its  weight-enumeator defined on T[n], i.e. the sum of y^(TotalHeight) over all trees in T[n]`):
print(``):
print(`Let R(x,y) be the bi-variate generating function: `):
print(``):
print(R(x,y)=Sum(P[n](y)*x^n,n=0..infinity)):
print(``):
print(`By generatingfunctionology,  R(x,y) satisfies the functional equation `):
print(``):
print(R(x,y)=x*(1+add( R(x*y,y)^i, i in DegSet))):
print(``):
print(`While R(x,y) does not seem to have closed-form, we can use it to compute the first`, N, `terms of the sequence of weight-enumerators`, P[n](y)):
print(``):
gu:=PolSeq(1+add(R^i,i in DegSet),R,y,N0):
print(gu):
print(``):
print(`and in Maple notation`):
print(``):
lprint(gu):
print(``):
print(`Let N[i](n) be the sum of the  i-th powers of the Total Heights of all trees in T[n]. Note that N[0](n) is the straight enumeration.`):
print(``):
print(`Let R[i](x) be the ordinary generating function of N[i](n), i.e. `):
print(``):	
print(R[i](x)=Sum(N[i](n)*x^n,n=0..infinity)):
print(``):	
print(`We have the following theorems expressing`, R[i](x), `in terms of R0(x) and x`):
print(``):
	
for i from 1 to k do
gu:=NiGF(1+add(R^i1, i1 in DegSet),R,i,x,R0):
print(`Theorem Number`, i, `: `):
print(``):
print(R[i](x)=subs(R0=R0(x),gu)):
print(``):
print(`and  in Maple notation, abbreviating R0(x) to R0 `):
print(``):
lprint(gu):
print(``):

od:


gu:=AlphaFFwNor(1+add(R1^i, i in DegSet),R1,N1,N2,k,MaxC):

print(`These can be used to crank-out many terms, obtain recurrences, that in turn generate MANY new terms, and using these`):
print(`we can estimate the limits of`):
print(``):
print(`the sequence of averages divided by n^(3/2)`):
print(``):
print(` from n=`, N1, `to n= `, N2, `are `):
print(``):
print(gu[1]):
print(``):

print(`the sequence of coefficients of variation (standard-deviation divided by the average) are `):
print(``):
print(`from n=`, N1, `to n= `, N2, `are `):
print(``):
print(gu[2]):
print(``):

for i from 3 to k do
print(`the sequence of coefficients of `, i, `scaled moments (aka alpha coefficients)  `):
print(``):
print(`from n= `, N1, `to n= `, N2, `are `):
print(``):
print(gu[i]):
print(``):
od:
print(``):
print(`----------------------------------------------------------------------------------`):
print(`This ends this article that took`, time()-t0, `seconds to generate. `):
print(`-------------------------------------------------------------------------------------`):

end:



#Paper2(DegSet,k,N0): Like Paper1(DegSet,k,MaxC,N0, N1,N2), but without the alpha coefficients
#For example, for rooted ordered trees where each vertex must have 0,1,or 2, children, type
#Paper2({1,2},4,20);
Paper2:=proc(DegSet,k,N0) local R,x,R0,gu,y,P,i,n,N,R1,t0,i1:

t0:=time():

print(``):
print(` Generating Functions for the Numerators of Moments of Total Height  on Ordered Trees where each vertex must have a Number of Children in the set`, DegSet union {0}):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`Let T(n) be the set of ordered (rooted) trees on n vertices where each vertex may be a leaf or else must have a number of children in the set`, DegSet):
print(``):
P:=1+add(R0(x)^i, i in DegSet):

print(`Let R0(x) be the ordinary generating function of the enumerating sequence. By generatingfunctionlogy, R0(x) is the unique formal power series`):
print(`satisfying `):
print(``):
print(R0(x)= x*P):
print(``):
print(`The first`, N0, `terms of the enumerating sequence starting at n=1, are `):

gu:=NumSeq(1+add(R1^i,i in DegSet),R1,N0):
print(``):
print(gu):
print(``):
print(`Let H[n] be the random variable "sum of heights (distances from the root) of all vertices", and let P[n](y) be the`):
print(`its  weight-enumeator defined on T[n], i.e. the sum of y^(TotalHeight) over all trees in T[n]`):
print(``):
print(`Let R(x,y) be the bi-variate generating function: `):
print(``):
print(R(x,y)=Sum(P[n](y)*x^n,n=0..infinity)):
print(``):
print(`By generatingfunctionology,  R(x,y) satisfies the functional equation `):
print(``):
print(R(x,y)=x*(1+add( R(x*y,y)^i, i in DegSet))):
print(``):
print(`While R(x,y) does not seem to have closed-form, we can use it to compute the first`, N, `terms of the sequence of weight-enumerators`, P[n](y)):
print(``):
gu:=PolSeq(1+add(R^i,i in DegSet),R,y,N0):
print(gu):
print(``):
print(`and in Maple notation`):
print(``):
lprint(gu):
print(``):
print(`Let N[i](n) be the sum of the  i-th powers of the Total Heights of all trees in T[n]. Note that N[0](n) is the straight enumeration.`):
print(``):
print(`Let R[i](x) be the ordinary generating function of N[i](n), i.e. `):
print(``):	
print(R[i](x)=Sum(N[i](n)*x^n,n=0..infinity)):
print(``):	
print(`We have the following theorems expressing`, R[i](x), `in terms of R0(x) and x`):
print(``):
	
for i from 1 to k do
gu:=NiGF(1+add(R^i1, i1 in DegSet),R,i,x,R0):
print(`Theorem Number`, i, `: `):
print(``):
print(R[i](x)=subs(R0=R0(x),gu)):
print(``):
print(`and  in Maple notation, abbreviating R0(x) to R0 `):
print(``):
lprint(gu):
print(``):

od:

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

end:



#NiGFbt(i,A): inputs a non-neg. integer i, and a symbol A, outputs an expression in A, let's call it P(A)
#such that NiGF(1+R^2,R,i,x,R0) equals P(A)/x, where A=sqrt(1-4*x^2):
#NiGFbt(2,A);
NiGFbt:=proc(i,A) local gu,R0,R,x:
gu:=x*NiGF(1+R^2,R,i,x,R0):
gu:=normal(subs(R0=(1-A)/2/x,gu)):
gu:=expand(subs(x=sqrt(1-A^2)/2,gu)):
gu:
end:


#NiBT(i,n): an explicit formula for the sum of the i-th powers of the total heights over all binary trees with 2n-1 vertices (and hence n leaves)
#here i is a specfic positive integer, and n is a symbol (denoting a postitive integer)
NiBT:=proc(i,n) local gu,gu1,gu2,A,r:

if i=0 then
 RETURN((2*n-2)!/n!/(n-1)!):
fi:

gu:=NiGFbt(i,A):
gu:=expand(subs(A=1/A,gu)):
gu1:=expand((gu+subs(A=-A,gu))/2):
gu2:=expand((gu-subs(A=-A,gu))/2):
gu1:=4^n*factor(add(coeff(gu1,A,2*r)*expand(binomial(r+n-1,r-1)),r=1..degree(gu1,A)/2)):
gu2:=factor(add(coeff(gu2,A,2*r+1)*simplify(convert(binomial(-1/2-r,n)/binomial(1/2,n),factorial)),r=0..trunc(degree(gu2,A)/2))):

gu1-2*((2*n-2)!/(n-1)!/n!)*gu2:

end:


#MomBT(i,n): Inputs a positive integer i, and outputs an explicit expression, in n, for the i-th moment about the mean
#(unless i=0 when it gives you the actual enumeration, and i=1 when it gives you the average) of the random variable
#"sum of distances from the root taken over all vertices" defined over all complete binary trees with n leaves (and hence 2n-1 vertices)
#Try:
#MomBT(2,n);
MomBT:=proc(i,n) local gu0,gu1,r:

gu0:=NiBT(0,n):
gu1:=NiBT(1,n):

if i=0 then
 RETURN(gu0):
fi:

if i=1 then
 RETURN(gu1/gu0):
fi:

add((-1)^r*binomial(i,r)*gu1^r*gu0^(i-r-1)*NiBT(i-r,n),r=0..i)/gu0^i:

end:


#NiBTc(i,n,cn): an explicit formula for the sum of the i-th powers of the total heights over all binary trees with 2n-1 vertices (and hence n leaves)
#divided by 4^n.
#Here i is a specfic positive integer, and n is a symbol (denoting a postitive integer) and cn is another symbol, short for (2n-2)!/(n-1)/n!
#Try:
#NiBTc(2,n,cn);
NiBTc:=proc(i,n,cn) local gu,gu1,gu2,A,r:

if i=0 then
 RETURN(cn):
fi:

gu:=NiGFbt(i,A):
gu:=expand(subs(A=1/A,gu)):
gu1:=expand((gu+subs(A=-A,gu))/2):
gu2:=expand((gu-subs(A=-A,gu))/2):
gu1:=4^n*factor(add(coeff(gu1,A,2*r)*expand(binomial(r+n-1,r-1)),r=1..degree(gu1,A)/2)):
gu2:=factor(add(coeff(gu2,A,2*r+1)*simplify(convert(binomial(-1/2-r,n)/binomial(1/2,n),factorial)),r=0..trunc(degree(gu2,A)/2))):

gu1-2*cn*gu2:


end:

#MomBTc(i,n,cn): Inputs a positive integer i, and a symbol cn, denoting (2n-2)!/(n!*(n-1)!)/4^n
#and outputs an explicit expression, in n, for the i-th moment about the mean
#(unless i=0 when it gives you the actual enumeration, and i=1 when it gives you the average) of the random variable
#"sum of distances from the root taken over all vertices" defined over all complete binary trees with n leaves (and hence 2n-1 vertices)
#Try:
#MomBTc(2,n,cn);
MomBTc:=proc(i,n,cn) local gu0,gu1,r:

gu0:=NiBTc(0,n,cn):
gu1:=NiBTc(1,n,cn):

if i=0 then
 RETURN(gu0):
fi:

if i=1 then
 RETURN(gu1/gu0):
fi:

add((-1)^r*binomial(i,r)*gu1^r*gu0^(i-r-1)*NiBTc(i-r,n,cn),r=0..i)/gu0^i:

end:

#AsyCn(n,K): the asymptotic expansion of binomial(2*n-2,n-1)/(4^n*n^(1/2)) up to order K. Try:
#AsyCn(n,5);
AsyCn:=proc(n,K) local lu:
 lu:=collect(expand(simplify(asympt((2*n-2)!/(n-1)!/n!/4^n/sqrt(n),n,K+1))),n);
 lu:=convert([op(1..nops(lu)-1,lu)],`+`):
 lu:
end:



#MomBTasy(i,m,K): Inputs a positive integer i, and outputs an ASYMPTOTIC expression, to order K, in m (that stands for sqrt(n)), for the i-th moment about the mean
#(unless i=0 when it gives you the actual enumeration, and i=1 when it gives you the average) of the random variable
#"sum of distances from the root taken over all vertices" defined over all complete binary trees with n leaves (and hence 2n-1 vertices)
#Try:
#MomBTasy(1,m,6);
MomBTasy:=proc(i,m,K) local lu,mu,i1,n:
lu:=collect(expand(simplify(asympt(MomBT(i,n)*sqrt(n),n,K))),n);

mu:=0:

for i1 from 1 to nops(lu) do
 if not type(op(i1,lu),function) then
  mu:=mu+op(i1,lu):
 fi:
od:

mu:=simplify(subs(n=m^2,mu),symbolic):
mu:=expand(mu/m):
mu:=collect(mu,m):

end:


#LimMomsBT(N): inputs a positive integer N, and outputs a list of numbers, of length N whose
#first component is the limit of of the average of the total height of complete binary trees with n leaves, divided by n^(3/2),
#The second component is the limit of the coeff. of variation and the third through the N componets are the
#limits of the scaled moments, aka as alpha coefficients. Try:
#LimMomsBT(4);
LimMomsBT:=proc(N) local gu,m,a,b,c,i:


a:=coeff(MomBTasy(1,m,5),m,3):
b:=coeff(MomBTasy(2,m,6),m,6):
gu:=[a,sqrt(b)/a]:

for i from 3 to N do

c:=coeff(MomBTasy(i,m,6+2*i),m,3*i):
gu:=[op(gu),c/b^(i/2)]:
od:

gu:

end:



#PaperBT(N): input a positive integer N and outputs an article about the explicit expressions, as well as asymptotic
#for the average, variance, and the third through the N-th moment, as well as the of the scaled limits
#for the random variable: total height defined over complete binary trees with n leaves. Try:
#PaperBT(4):
PaperBT:=proc(N) local n,m,gu,t0,i:

t0:=time():

print(``):
print(`On the Statistical Distribution of the Total Height  on Complete (Ordered) Binary Trees `):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`Let B(n) be the set of ordered (rooted) trees on n  LEAVES where each vertex may be a leaf or else must have exactly two children`):

gu:=MomBT(1,n):

print(`Theorem 1: The average is given explicity by`):
print(``):
print(gu):
print(``):
print(`and in Maple notation by`):
print(``):
lprint(gu):
print(``):
gu:=MomBTasy(1,m,10):
print(`this is asymptotically, where m=sqrt(n) `, gu):
print(``):
print(`and in Maple notation`):
print(``):
lprint(gu):
print(``):

gu:=MomBT(2,n):
print(`Theorem 2: The variance is given explicity by`):
print(``):
print(gu):
print(``):
print(`and in Maple notation by`):
print(``):
lprint(gu):
print(``):
gu:=MomBTasy(2,m,10):
print(`this is asymptotically, where m=sqrt(n) `, gu):
print(``):
print(`and in Maple notation`):
print(``):
lprint(gu):
print(``):

for i from 3 to N do

gu:=MomBT(i,n):
print(`Theorem Number`, i, `: The `, i, `-th moment (about the mean) is given explicity by`):
print(``):
print(gu):
print(``):
print(`and in Maple notation by`):
print(``):
lprint(gu):
print(``):
gu:=MomBTasy(i,m,10+i):
print(`this is asymptotically, where m=sqrt(n) `, gu):
print(``):
print(`and in Maple notation`):
print(``):
lprint(gu):
print(``):
od:

gu:=LimMomsBT(N):
print(`Finally, the limit of the average divided by n^(3/2), the limit of the coefficient of variation, and the limit of the scaled moment up to the`, N, `-th `):
print(`are ` ):
print(``):
print(gu):
print(``):
print(`and in Maple notation`):
lprint(gu):
print(`And in floating-point: `):
print(``):
lprint(evalf(evalf(gu,100),20)):
print(``):

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