%I #16 Nov 02 2019 08:45:40
%S 1,1,1,2,8,46,337,2976,30627,359222,4725968,68903766,1102712316,
%T 19219507328,362428546833,7352854216056,159705991698432,
%U 3697928742242694,90933523698184947,2366758931071064064
%N Coefficient of x^n in A(x)^n is A000670(n), which gives preferential arrangements of n labeled elements.
%H Vaclav Kotesovec, <a href="/A088791/b088791.txt">Table of n, a(n) for n = 0..333</a>
%F G.f. satisfies: A(x)^2 = A(x*A(x)) + x*A(x).
%F a(n) ~ (n-1)! / (4 * (log(2))^(n+1)). - _Vaclav Kotesovec_, Feb 12 2015
%F O.g.f.: A(x) = x/( series reversion x*B(x) ), where B(x) = 1 + x + 2*x^2 + 6*x^3 + 25*x^4 + 137*x^5 + ... is the o.g.f. of A084784. - _Peter Bala_, Jun 23 2015
%t nmax = 19; sol = {a[0] -> 1};
%t Do[A[x_] = Sum[a[k] x^k, {k, 0, n}] /. sol; eq = CoefficientList[A[x]^2 - (A[x A[x]] + x A[x]) + O[x]^(n + 1), x] == 0 /. sol; sol = sol ~Join~ Solve[eq][[1]], {n, 1, nmax}];
%t sol /. Rule -> Set;
%t a /@ Range[0, nmax] (* _Jean-François Alcover_, Nov 02 2019 *)
%o (PARI) {a(n)=local(A, m); if(n<1, n==0, m=1; A=1+x; for(i=1, n, A=(subst(A, x, x*A+x*O(x^n)) + x*A)/A); polcoeff(A, n))}
%o for(n=0, 20, print1(a(n), ", ")) \\ _Vaclav Kotesovec_, Feb 11 2015
%Y Cf. A000670, A084784.
%K nonn,easy
%O 0,4
%A _Paul D. Hanna_, Oct 16 2003