%I #9 Aug 14 2013 15:31:49
%S 1,1,3,13,97,861,10171,144313,2425473,46361017,1008845011,24440301381,
%T 653993215393,19126571703253,607566772915467,20816075734498801,
%U 765497764431847681,30064774690536609393,1256227494273614356003,55637289075570248646397,2603702479493046357670881
%N Number of sets (forests) of labeled identity trees (trees enumerated by A228159) with n nodes.
%H Alois P. Heinz, <a href="/A228160/b228160.txt">Table of n, a(n) for n = 0..150</a>
%F E.g.f.: exp(A(x)) where A(x) is the e.g.f. for A228159.
%p b:= proc(n, i) option remember; `if`(n=0, 1, `if`(i<1, 0,
%p add(binomial(b(i-1$2), j)*b(n-i*j, i-1), j=0..n/i)))
%p end:
%p a:= proc(n) option remember; `if`(n=0, 1,
%p add(binomial(n-1, j-1) *j!*b(j-1$2)*a(n-j), j=1..n))
%p end:
%p seq(a(n), n=0..25); # _Alois P. Heinz_, Aug 14 2013
%t nn = 20; SolveAlways[
%t 0 == Series[
%t f[x] - x Product[(1 + x^i)^(a[i]/i!), {i, 1, nn}], {x, 0,
%t nn}], x]; b = Sum[a[n] x^n/n!, {n, 1, nn}] /. sol;
%t Range[0, nn]! Flatten[CoefficientList[Series[Exp[b], {x, 0, nn}], x]]
%K nonn
%O 0,3
%A _Geoffrey Critzer_, Aug 14 2013