%I #17 Sep 08 2022 08:46:05
%S 0,0,0,1,16,177,1726,15912,143148,1279939,11504326,104686659,
%T 968808308,9144180028,88184565504,869867691833,8781919559956,
%U 90765497635245,960434143555986,10403548856756708,115336464546432180,1308260884070774299,15177980646442995698,180036437138753006607,2182526416321158803528
%N 4*B(n+4) - (4*n+15)*B(n+3) + (n^2+8*n+9)*B(n+2) - (4*n+3)*B(n+1) + n*B(n), where B(i) are the Bell numbers A000110.
%H Vincenzo Librandi, <a href="/A226507/b226507.txt">Table of n, a(n) for n = 0..200</a>
%H B. Chern, P. Diaconis, D. M. Kane, and R. C. Rhoades, <a href="https://arxiv.org/abs/1304.4309">Closed expressions for averages of set partition statistics</a>, arXiv:1304.4309 [math.CO], 2013.
%F a(n) ~ n^4 * Bell(n) / LambertW(n)^2 * (1 - 4/LambertW(n) + 4/LambertW(n)^2). - _Vaclav Kotesovec_, Jul 28 2021
%p A000110 := proc(n) option remember; if n <= 1 then 1 else add( binomial(n-1, i)*A000110(n-1-i), i=0..n-1); fi; end;
%p B:=A000110;
%p f:=n->4*B(n+4) - (4*n+15)*B(n+3) + (n^2+8*n+9)*B(n+2) - (4*n+3)*B(n+1) + n*B(n);
%p [seq(f(n),n=0..30)];
%t Table[4 BellB[n+4] - (4 n + 15) BellB[n + 3] + (n^2 + 8 n + 9) BellB[n+2] - (4 n + 3) BellB[n+1] + n BellB[n],{n, 0, 30}] (* _Vincenzo Librandi_, Jul 16 2013 *)
%o (Magma) [4*Bell(n+4)-(4*n+15)*Bell(n+3)+(n^2+8*n+9)*Bell(n+2)-(4*n+3)*Bell(n+1)+n*Bell(n): n in [0..30]]; // _Vincenzo Librandi_, Jul 16 2013
%Y Cf. A000110, A200580.
%K nonn
%O 0,5
%A _N. J. A. Sloane_, Jun 10 2013