# Maple program for enumeration of rooted trees by height
# (after Riordan's IBM J. paper)

did:=proc(m,n) if irem(m,n) = 0 then 1           else 0: fi: end:
interface(screenwidth=1000);

# EULER takes [a_1,a_2,a_3,...] and produces [b_1,b_2,b_3,...] with
# 1 + Sum_{n=1..inf} b_n x^n = 1 / Product_{n=1..inf} (1-x^n)^a_n.
EULER:=proc(a)  local b,c,i,d:
if whattype(a) <> list then RETURN([]); fi:
c:=[]:
for i to nops(a) do c:=[op(c), add( d*did(i,d)*a[d] , d=1..i)]: od:
b:=[]:
for i to nops(a) do
b:=[op(b),(1/i)*(c[i]+ add( c[d]*b[i-d], d=1..i-1))]: od:
RETURN(b);
end:

M:=102; Mh:=102; Mp:=20;

R[0]:=x;
R[1]:=x^2/(1-x);

r[1,1]:=0: for n from 2 to M do r[1,n]:=1; od:
for n from 1 to M do s[1,n]:=1; od:

for h from 2 to Mh do
a:=[seq(s[h-1,n],n=1..M)];
b:=[1, op(EULER(a))];
for n from 1 to M do s[h,n]:=b[n]; od:
for n from 1 to M do r[h,n]:=s[h,n] - s[h-1,n]; od:
print("b file for r = trees of height ", h);
print([seq(r[h,n],n=1..Mp)]);
#for n from 1 to M do lprint(n, r[h,n]); od:
print("b file for s = trees of height <= ", h);
print([seq(s[h,n],n=1..Mp)]);
#for n from 1 to M do lprint(n, s[h,n]); od:
R[h]:=add(r[h,n]*x^n,n=1..M);
                   od:

t1:=add( (1/2)*(R[h]^2+subs(x=x^2,R[h])),h=0..Mh):
t2:=series(t1,x,2*Mh):
for n from 0 to 2*Mh-1 do lprint(n, coeff(t2,x,n)); od: