######################################################################
## Krandick.txt Save this file as  Krandick.txt  to use it,          #
# stay in the                                                        #
## same directory, get into Maple (by typing: maple <Enter> )        #
## and then type:  read  `Krandick.txt` <Enter>                      #
## Then follow the instructions given there                          #
##                                                                   #
## Written by Doron Zeilberger, Rutgers University ,                 #
##  DoronZeil at gmail dot com                                       # 
######################################################################


print(`First Written: April  2024: tested for Maple 2020 `):
print(`Version  : April 26, 2024 `):
print():
print(`This is Krandick.txt, A Maple package`):
print(`accompanying Shalosh B. Ekhad and Doron Zeilberger's  article: `):
print(` Exploring Werner Krandick's Binary Tree Statistics`):
print(` http://sites.math.rutgers.edu/~zeilberg/tokhniot/mamarim/mamarimhtml/krandick.html `):
print():
print(`The most current version is available on WWW at:`):
print(` http://sites.math.rutgers.edu/~zeilberg/tokhniot/Krandick.txt .`):
print(`Please report all bugs to: DoronZeil at gmail dot com .`):
print():
print(`For general help, and a list of the MAIN functions,`):
print(` type "ezra();". For specific help type "ezra(procedure_name);" `):
print(`For a list of the supporting functions type: ezra1();`):
print():
print(`-------------------`):

print(`For a list of the generating functions procedures type: ezraGF();`):

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





ezraGF:=proc()
if args=NULL then

print(`The  generating function  procedures are: `):
print(` Bx, Dxq, Jxq, Jxtq,  ProveJxtq, Rxt `):

else
ezra(args):

fi:


end:

ezraS:=proc()
if args=NULL then

print(`The Story procedures are: `):
print(``):

else
ezra(args):

fi:


end:

ezra1:=proc()
if args=NULL then

print(`The SUPPORTING procedures are:  Alpha, DFSt, Empir, GR, LD, NJ, NJD, NuIV, RA, SeqJ,`):
print(``):

else
ezra(args):

fi:


end:

ezra:=proc()
if args=NULL then

print(` Krandick.txt: A Maple package for  studying Wener Krandick's binary tree statistics`):

print(`The MAIN procedures are: BT, BTd,  Dnq,  Jnq, Jntq,RBT,  Rnt `):


print(``):


elif nargs=1 and args[1]=Alpha then
print(`Alpha(f,x,N): given a p.g.f. outputs the average, variance, and scaled moments up to the N-th. Try`):
print(`Alpha((1+x)^100,x,4);`):


elif nargs=1 and args[1]=BT then
print(`BT(n): the set of full binary trees with n internal vertices (and hence altogether 2n+1 vertices) labeled in [] notation. Try:`):
print(`BT(4);`):

elif nargs=1 and args[1]=BTd then
print(`BTd(n): the set of full binary trees with n internal vertices (and hence altogether 2n+1 vertices) labeled according to Depth First Search. Try:`):
print(`BTd(4);`):

elif nargs=1 and args[1]=Bx then
print(`Bx(x): The generating function whose coefficient of x^n is the number of full binary trees with n internal nodes, i.e. the famous Cataln numbers.`):
print(`Try:Bx(x);`):


elif nargs=1 and args[1]=Dxq then
print(`Dxq(x,q): the generating function whose coefficient of x^n*q^i is the number of full binary trees with n internal nodes and total sum distance i. Try:`):
print(`Dxq(x,q);`):


elif nargs=1 and args[1]=Dnq then
print(`Dnq(n,q): the weight-enumerator of full binary trees with n internal leaves according to the weight q^TotalumpDistances. Try`):
print(`Dnq(10,q);`):

elif nargs=1 and args[1]=DFSt then
print(`DFSt(T): inputs a full binary tree T, in [] notation,  with n internal vertices, say, and outputs it as a list, where the vertices are labelled {1,...2*n+1}  and L[i] is []`):
print(`if it is a leaf or [a1,a2] if a1 and a2 are its left and right children. Try:`):
print(`DFSt([[],[]]):`):


elif nargs=1 and args[1]=Empir then
print(`Empir(L,x,P): inputs a  sequence of numbers L and tries to guess an equation of the form F(x,P)=0, where F is a bi-variate polynomial satisfied by`):
print(`Sum(L[i+1]*x^i,i=0..infinity). Try:` );
print(`Empir([seq(binomial(2*i,i)/(i+1),i=0..20)],x,P);`):


elif nargs=1 and args[1]=Jxq then
print(`Jxq(x,q): the generating function whose coefficient of x^n*q^i is the number of full binary trees with n internal nodes and total sum distance i. Try:`):
print(`Jxq(x,q);`):

elif nargs=1 and args[1]=Jxtq then
print(`Jxtq: the generating function whose coefficient of x^n*q^i*t^j is the number of full binary trees with n internal nodes i and number of jumps i and distance of rightmost root j. Try:`):
print(`Jxtq(x,t,q);`):

elif nargs=1 and args[1]=Jnq then
print(`Jnq(n,q): the weight-enumerator of full binary trees with n internal leaves according to the weight q^Number of Jumps. Try`):
print(`Jnq(10,q);`):



elif nargs=1 and args[1]=Jntq then
print(`Jntq(n,t,q): the weight-enumerator of full binary trees with n internal leaves according to the weight t^DepthOfRightMostLeaf*q^Number of Jumps. Try: `):
print(`Jntq(5,t,q);`):

elif nargs=1 and args[1]=GR then
print(`GR(L,n): Inputs a list of numbers L and a variable n and outputs the rational function  of degree(numer)+degree(denom) <=nops(L)-2`):
print(`RP(n), such that R(i)=L[i], for all i from 1 to nops(L). OR FAIL`):

elif nargs=1 and args[1]=LD then
print(`LD(L): Given a list of integers L outputs i with prob. L[i]/(L[1]+..+L[n]). Try:`):
print(`LD([1,2,2,1]);`):

elif nargs=1 and args[1]=NJ then
print(`NJ(T): the number of jumps of the full binary tree T`):


elif nargs=1 and args[1]=NJD then
print(`NJD(T): the total of jump distances of the full binary tree T`):

elif nargs=1 and args[1]=NuIV then
print(`NuIV(T): the number of internal vertices of the full binary tree T. Try:`):
print(`NuIV([[],[]]):`):

elif nargs=1 and args[1]=ProveJxtq then
print(`ProveJxtq():  proves the conjectured expression for Jxtq(x,t,q). Try:`):
print(`ProveJxtq();`):

elif nargs=1 and args[1]=RA then
print(`RA(T): the length of the rightmost branch of the full binary tree T. Try:`):
print(`RA([[],[]]);`):

elif nargs=1 and args[1]=RBT then
print(`RBT(n): a random binary tree with n internal vertices. Try:`):
print(`RBT(5);`):


elif nargs=1 and args[1]=Rnt then
print(`Rnt(n,t): the weight-enumerator of full binary trees according to the t to the depth of the rightmost leaf. Try:`):
print(`Rnt(4,t);`):

elif nargs=1 and args[1]=Rxt then
print(`Rxt(x,t): The generating function whose coefficient of x^n*t^i is the number of full binary trees with n internal vertices and distance i from the root to the rightmost tree`):
print(`in other words Sum(Rnt(n,t)*x^n,n=0..infinity). Try:`):
print(`Rxt(x,t);`):


elif nargs=1 and args[1]=SeqJ then
print(`SeqJ(t,q,N): same as [seq(Jntq(n,t,q),n=0..N)] but via the functional equation. Try:`):
print(`SeqJ(t,q,10);`):


else
 print(`There is no such thing as`, args):

fi:


end:



#start from SCHUTZENBERGER

#empir(gu,degx,degP,x,P) 
#to "fit" an algebraic equation F(P(x),x)=0 of degree
#degP in P(x) and degx in n for P(x):=sum_i gu[i]*x^i
 
empir:=proc(gu,degx,degP,x,P)
local i1,i2,F,a,cand,lu,eq,var,mu,flo,pip:
if (1+degx)*(1+degP) > nops(gu)-3 then
 RETURN(`sequence too small`):
fi:
 
F:=0:
var:={}:
for i1 from 0 to degx do
 for i2 from 0 to degP do
  F:=F+a[i1,i2]*x^i1*P^i2:
  var:=var union {a[i1,i2]}:
 od:
od:
 
cand:=0:
 
for i1 from 0 to nops(gu)-1 do
 cand:=cand+op(i1+1,gu)*x^i1:
od:
 
lu:=subs(P=cand,F):
lu:=taylor(lu,x=0,nops(gu)-1):
 
eq:={}:
 
for i1 from 0 to nops(gu)-2 do
eq:=eq union {coeff(lu,x,i1)=0}
od:
 
mu:=solve(eq,var):
 
 
F:=subs(mu,F):
 
if F=0 then
 RETURN(0):
fi:
 
flo:=degree(F,P):
pip:=coeff(F,P,flo):
flo:=degree(pip,x):
pip:=coeff(pip,x,flo):
F:=F/pip:
normal(F):
end:

#Empir(gu,x,P)
Empir:=proc(gu,x,P)
local degx,degP,L,lu:
 
for L from 1 to (nops(gu)-3)/3 do
for degP from 1 to L do
for degx from 0 to min(trunc((nops(gu)-3)/(1+degP))-1,L-degP) do
 
lu:=empir(gu,degx,degP,x,P):
 
if lu<>0 then
RETURN(lu):
fi:
od:
od:
od:
0:
end:
 
#end from SCHUTZENBERGER

#GR1(L,n,d1,d2): Inputs a list of numbers L and a variable n and a non-neg. d outputs the rational function of degree d
#P(n), such that P(i)=L[i], for all i from 1 to nops(L). OR FAIL
#nops(L)>=d+2
GR1:=proc(L,n,d1,d2) local a,i,P,eq,var,b:
if nops(L)<d1+d2+4 then
 ERROR(`List too short`):
 fi:

P:=add(a[i]*n^i,i=0..d1)/(n^d2+add(b[i]*n^i,i=0..d2-1)):


var:={seq(a[i],i=0..d1),seq(b[i],i=0..d2-1)}:
eq:={seq(numer(L[i]-subs(n=i,P)),i=1..d1+d2+4)}:
var:=solve(eq,var):

 if var=NULL then
   RETURN(FAIL):
 fi:
P:=subs(var,P):
if {seq(L[i]-subs(n=i,P),i=d1+d2+4..nops(L))}<>{0} then
  RETURN(FAIL):
 fi:
normal(P):
end:



#GR(L,n): Inputs a list of numbers L and a variable n and outputs the polynomial of degree <=nops(L)-2
#P(n), such that P(i)=L[i], for all i from 1 to nops(L). OR FAIL
#nops(L)>=d+2
GR:=proc(L,n) local d,d1,P:
for d from 1 to (nops(L)-4)/2 do
for d1 from 1 to d do
 P:=GR1(L,n,d1,d-d1):
 if P<>FAIL then
    RETURN(P):
 fi:
od:
od:

FAIL:
end:


##Start from HISTABRUT
#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:


#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(gu):
fi:

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

end:


##end from HISTABRUT

#BTd(n): the set of full binary trees with n internal vertices (and hence altogether 2n+1 vertices) labeled according to Depth First Search. Try:
#BTd(3);
BTd:=proc(n) local gu,gu1,gu2,k,i,gu11,gu21:
option remember:
if n=0 then
 RETURN({{}}):
fi:
gu:={}:

for k from 0 to n-1 do
gu1:=BTd(k):
gu2:=BTd(n-1-k):
gu1:=subs({seq(i=i+1,i=1..3*n)},gu1):
gu2:=subs({seq(i=i+2*k+2,i=1..3*n)},gu2):
gu:=gu union {seq(seq({{1,2},{1,2*k+3},op(gu11),op(gu21)},gu11 in gu1),gu21 in gu2)}:
od:
gu:
end:


#Rnt(n,t): The weight enumerator of full binary trees with n internal vertices according to the depth of the rightmost vertex. Try:
#Rnt(n,t)
Rnt:=proc(n,t) local k:
option remember:
if n=0 then
 RETURN(1):
fi:
expand(add(t*subs(t=1,Rnt(k,t))*Rnt(n-1-k,t),k=0..n-1)):
end:

#Jntq(n,t,q): the weight-enumerator of full binary trees with n internal leaves according to the weight t^DepthOfRightMostLeaf*q^Number of Jumps
Jntq:=proc(n,t,q) local k,gu,lu,i:
option remember:

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

gu:=0:

for k from 0 to n-1 do
 lu:=Jntq(k,t,q):
gu:=expand(gu+t*coeff(lu,t,0)*Jntq(n-1-k,t,q)+ q*t*add(coeff(lu,t,i),i=1..degree(lu,t))*Jntq(n-1-k,t,q)):
od:
gu:
end:


#Jnq(n,q): the weight-enumerator of full binary trees with n internal leaves according to the weight q^NumberOfJumps
Jnq:=proc(n,q) local t:subs(t=1,Jntq(n,t,q)):end:


#Dnq(n,q): the weight-enumerator of full binary trees with n internal leaves according to the weight q^TotalJumpDistances
Dnq:=proc(n,q) local t: expand(q^n*subs(t=1/q,Rnt(n,t))):end:





#SeqJ(t,q,N): same as [seq(Jntq(n,t,q),n=0..N)] but via the functional equation. Try:
#SeqJ(t,q,10);
SeqJ:=proc(t,q,N) local phi,phi1,i,x:
phi:=1:
phi1:=expand(1+x*t*subs(t=0,phi)*phi+x*t*q*(subs(t=1,phi)-subs(t=0,phi))*phi):
phi1:=add(coeff(phi1,x,i)*x^i,i=0..N):

while phi<>phi1 do
 phi:=phi1:
 phi1:=expand(1+x*t*subs(t=0,phi)*phi+x*t*q*(subs(t=1,phi)-subs(t=0,phi))*phi):
 phi1:=add(coeff(phi1,x,i)*x^i,i=0..N):
od:
[seq(coeff(phi,x,i),i=0..N)]:
end:


#SeqD(t,q,N): same as [seq(Dntq(n,t,q),n=0..N)] but via the functional equation. Try:
#SeqD(t,q,10);
SeqD:=proc(t,q,N) local phi,phi1,i,x:
phi:=1:
phi1:=expand(1+x*t*subs(t=q,phi)*phi):
phi1:=add(coeff(phi1,x,i)*x^i,i=0..N):

while phi<>phi1 do
 phi:=phi1:
phi1:=expand(1+x*t*subs(t=q,phi)*phi):
 phi1:=add(coeff(phi1,x,i)*x^i,i=0..N):
od:
[seq(coeff(phi,x,i),i=0..N)]:
end:



#BT(n): the set of full binary trees in [] notation. Try
#BT(3);
BT:=proc(n) local k,gu1,gu2,T1,T2,gu :
option remember:

if n=0 then
 RETURN({[]}):
fi:

gu:={}:

for k from 0 to n-1 do
 gu1:=BT(k):
 gu2:=BT(n-1-k):
 gu:= gu union {seq(seq([T1,T2],T1 in gu1), T2 in gu2)}:
od:
gu:
end:





#LD(L): Given a list of integers L outputs i with prob. L[i]/(L[1]+..+L[n]). Try:
#LD([1,2,2,1]);
LD:=proc(L) local n,i,M,T,c,k:
n:=nops(L):
T:=convert(L,`+`):
M:=[]:

c:=0:
for i from 1 to n do
 c:=c+L[i]:
 M:=[op(M),c]:
od:

k:=rand(1..T)():

for i from 1 while k>M[i] do od:
i:

end:

 

#RBT(n): a random binary tree with n internal vertices. Try:
#RBT(5);
RBT:=proc(n) local k,L,T1,T2:
if n=0 then
 RETURN([]):
fi:
L:=[seq(binomial(2*k,k)/(k+1)*binomial(2*(n-1-k),(n-1-k))/(n-k),k=0..n-1)]:

k:=LD(L)-1:

T1:=RBT(k):
T2:=RBT(n-1-k):
[T1,T2]:


end:

#RA(T): the length of the rightmost branch of the full binary tree T. Try:
#RA([[],[]]);
RA:=proc(T):if T=[] then 0: else 1+RA(T[2]):fi: end:

#NJ(T): the number of jumps of the full binary tree T
NJ:=proc(T):

if T=[] then
 RETURN(0):
fi:

if T[1]=[] then
 RETURN (NJ(T[2])):
else
RETURN (NJ(T[1])+NJ(T[2])+1):
fi:

end:


#NJD(T): the number of total jump distances of the full binary tree T
NJD:=proc(T):
NuIV(T)-RA(T):
end:

#NuIV(T): the number of internal vertices of the full binary tree T. Try:
#NuIV([[],[]]):
NuIV:=proc(T):
if T=[] then
 0:
else
NuIV(T[1])+NuIV(T[2])+1:
fi:
end:



#DFSt(T): inputs a full binary tree T, in [] notation,  with n internal vertices, say, and outputs it as a list, where the vertices are labelled {1,...2*n+1}  and L[i] is []
#if it is a leaf or [a1,a2] if a1 and a2 are its left and right children. Try:
#DFSt([[],[]]):
DFSt:=proc(T) local k1,k2,L1,L2,i:
if T=[] then
 RETURN([[]]):
fi:
k1:=NuIV(T[1]):
k2:=NuIV(T[2]):



L1:=DFSt(T[1]):
L1:=subs({seq(i=i+1,i=1..2*k1+1)},L1):

L2:=DFSt(T[2]):

L2:=subs({seq(i=i+2*k1+2,i=1..2*k2+1)},L2):

[[2,2*k1+3],op(L1),op(L2)]:
end:




###start generating function section

#Bx(x): The generating function whose coefficient of x^n is the number of full binary trees with n internal nodes, i.e. the famous Cataln numbers.
#Try:
#Bx(x);
Bx:=proc(x): (1-sqrt(1-4*x))/(2*x):end:

#Rxt(x,t): The generating function whose coefficient of x^n*t^i is the number of full binary trees with n internal vertices and distance i from the root to the rightmost tree
#in other words Sum(Rnt(n,t)*x^n,n=0..infinity). Try:
#Rxt(x,t);
Rxt:=proc(x,t):
normal(1/(1-x*t*Bx(x))):
end:

#Dxq(x,q): the generating function whose coefficient of x^n*q^i is the number of full binary trees with n internal nodes and total sum distance i. Try:
#Dxq(x,q);
Dxq:=proc(x,q) local t:
normal(subs({x=q*x,t=1/q},Rxt(x,t))):
end:

#Jxq: the generating function whose coefficient of x^n*q^i is the number of full binary trees with n internal nodes i and number of jumps i. Try:
#Jxq(x,q);
Jxq:=proc(x,q): 
-1/2*(-q*x+(q^2*x^2-2*q*x^2-2*q*x+x^2-2*x+1)^(1/2)+x-1)/q/x:
end:

#Jxtq: the generating function whose coefficient of x^n*q^i*t^j is the number of full binary trees with n internal nodes i and number of jumps i and distance of rightmost root j. Try:
#Jxtq(x,q);
Jxtq:=proc(x,t,q): 
-1/2*(-q*t*x+t*x+(q^2*t^2*x^2-2*q*t^2*x^2-2*q*t^2*x+t^2*x^2-2*t^2*x+t^2)^(1/2)+t-2)/(q*t*x+t^2*x-t*x-t+1):
end:

#ProveJxtq():  proves the conjectured expression for Jxtq(x,t,q). Try:
#ProveJxtq();
ProveJxtq:=proc() local f,x,t,q,gu: 
f:=Jxtq(x,t,q): 
gu:=f-1-x*t*subs(t=0,f)*f-x*t*q*(subs(t=1,f)-subs(t=0,f))*f:
evalb(simplify(gu,symbolic)=0):
end:

###end generating function section