%I #27 Sep 08 2022 08:46:05
%S 0,16,160,1686,21276,328498,6149136,137105016,3577543452,107601726030,
%T 3683660206080,142035221781402,6113719409724768,291540411275223912,
%U 15300594717301253800,878667035554110785662,54932693182800769213284
%N a(n) = B(n+3) - 6*B(n+2) + 8*B(n+1)*B(n+1) - B(n), where the B(i) are Bell numbers (A000110).
%H Vincenzo Librandi, <a href="/A225591/b225591.txt">Table of n, a(n) for n = 0..200</a>
%H B. Chern, P. Diaconis, D. M. Kane, R. C. Rhoades, <a href="http://math.stanford.edu/~rhoades/FILES/setpartitions.pdf">Closed expressions for averages of set partition statistics</a>, 2013.
%t Table[BellB[n+3] - 6 BellB[n+2] + 8 BellB[n+1] BellB[n+1] - BellB[n], {n, 0, 20}] (* _Vincenzo Librandi_, Jul 16 2013 *)
%t #[[4]]-6#[[3]]+8#[[2]]^2-#[[1]]&/@Partition[BellB[Range[0,20]],4,1] (* _Harvey P. Dale_, Nov 01 2016 *)
%o (PARI) B(n) = if (n<=1, return (1), return (sum(i=0, n-1, binomial(n-1, i)*B(n-1-i))))
%o a(n) = B(n+3) - 6*B(n+2) + 8*B(n+1)*B(n+1) - B(n)
%o (Magma) [Bell(n+3)-6*Bell(n+2)+8*Bell(n+1)*Bell(n+1)-Bell(n): n in [0..20]]; // _Vincenzo Librandi_, Jul 16 2013
%Y Cf. A005493, A226506 (see Prop 3.1 (i) in Chern et al. link).
%K nonn
%O 0,2
%A _Michel Marcus_, Jun 19 2013