OFFSET
1,6
COMMENTS
From R. J. Mathar, Apr 13 2019: (Start)
The associated triangle H_{p,j}, p >= 1, 1 <= j <= p, a(n) = Sum_{j=1..p} j*H_{p,j}, row sums in A001678, starts:
1;
0, 0;
1, 0, 0;
1, 0, 0, 0;
1, 0, 1, 0, 0;
1, 1, 1, 0, 0, 0;
2, 2, 1, 0, 1, 0, 0;
1, 4, 2, 2, 1, 0, 0, 0;
3, 4, 4, 5, 2, 0, 1, 0, 0;
3, 7, 7, 9, 4, 4, 1, 0, 0, 0;
5, 9, 15, 14, 11, 9, 3, 0, 1, 0, 0;
4, 14, 23, 28, 25, 19, 7, 6, 1, 0, 0, 0;
11, 15, 39, 46, 55, 38, 24, 14, 5, 0, 1, 0, 0;
6, 32, 54, 86, 97, 86, 64, 36, 11, 9, 1, 0, 0, 0;
(End)
REFERENCES
N. J. A. Sloane and Simon Plouffe, The Encyclopedia of Integer Sequences, Academic Press, 1995 (includes this sequence).
LINKS
F. Harary and E. M. Palmer, Probability that a point of a tree is fixed, Math. Proc. Camb. Phil. Soc. 85 (1979) 407-415.
MAPLE
Hpj := proc(Hofxy, p, j)
coeftayl(Hofxy, x=0, p) ;
coeftayl(%, y=0, j) ;
simplify(%) ;
end proc:
Hxy := proc(x, y, pmax, hxyinit)
if pmax = 0 then
x*y ;
else
pp := 1;
for p from 1 to pmax do
t :=1 ;
for j from 1 to p do
t := t*(1+x^p*y^j+add(x^(k*p), k=2..pmax+1))^Hpj(hxyinit, p, j) ;
end do:
pp := pp*t ;
end do:
x*y*%/(1+x*y) ;
end if;
end proc:
hxy := Hxy(x, y, 0, 0) ;
for pmax from 2 to 20 do
Hxy(x, y, pmax, hxy) ;
taylor(%, x=0, pmax+2) ;
convert(%, polynom) ;
taylor(%, y=0, pmax+2) ;
hxy := convert(%, polynom) ;
for p from 0 to pmax do
Ap := 0 ;
for j from 1 to p do
Ap := Ap+j*Hpj(hxy, p, j) ;
end do:
printf("%d, ", Ap) ;
end do:
print() ;
end do: # R. J. Mathar, Apr 13 2019
MATHEMATICA
Hpj[Hofxy_, p_, j_] := SeriesCoefficient[SeriesCoefficient[Hofxy, {x, 0, p}] , {y, 0, j}];
Hxy [x_, y_, pMax_, hxyinit_] := If [pMax == 0, x y, pp = 1; For[p = 1, p <= pMax, p++, t = 1; For[j = 1, j <= p, j++, t = t(1 + x^p y^j + Sum[x^(k*p), {k, 2, pMax + 1}])^Hpj[hxyinit, p, j]]; pp = pp t]; x*y* pp/(1 + x y)];
hxy = Hxy[x, y, 0, 0];
Reap[For[pMax = 2, pMax <= terms - 1, pMax++, Print["pMax = ", pMax]; sx = Series[Hxy[x, y, pMax, hxy], {x, 0, pMax + 2}] // Normal; sy = Series[sx, {y, 0, pMax + 2}]; hxy = sy // Normal; For[p = 0, p <= pMax, p++, Ap = 0; For[j = 1, j <= p, j++, Ap = Ap + j Hpj[hxy, p, j]]; If[pMax == terms - 1, Print[Ap]; Sow[Ap]]]]][[2, 1]] (* Jean-François Alcover, Mar 22 2020, after R. J. Mathar *)
CROSSREFS
KEYWORD
nonn,easy
AUTHOR
EXTENSIONS
More terms from R. J. Mathar, Apr 13 2019
STATUS
approved