%I #46 Dec 14 2023 17:26:11
%S 1,1,6,48,408,5040,72000,1184400,22619520,482993280,11459750400,
%T 299495750400,8531976499200,263353163673600,8754879893760000,
%U 311808414677760000,11845876873678848000,478163414336864256000,20436460099541950464000,921972301728418676736000
%N Expansion of e.g.f. 1/(1-x-2*x^2-3*x^3).
%C a(n) is the number of ways to partition [n] into blocks of size at most 3, order the blocks, order the elements within each block, and choose 1 element from each block.
%C E.g.: a(4) = 408 since we have the following cases:
%C 1,2,3,4: 24 such orderings, 1 way to choose one element from each block;
%C 12,34: 24 such orderings, 2*2 ways to choose one element from each block;
%C 12,3,4: 72 such orderings, 2*1*1 ways to choose one element from each block;
%C 123,4: 48 such orderings, 3*1 ways to choose one element from each block;
%C so 24*1 + 24*4 + 72*2 + 48*3 = 408 ways.
%F a(n) = A000142(n)*A101822(n).
%F a(n) = n*(a(n-1)+(n-1)*(2*a(n-2)+(n-2)*3*a(n-3))) for n>=3. - _Alois P. Heinz_, Dec 14 2023
%p a:= proc(n) option remember; `if`(n=0, 1, add(
%p a(n-j)*binomial(n, j)*j!*j, j=1..min(3, n)))
%p end:
%p seq(a(n), n=0..19); # _Alois P. Heinz_, Dec 14 2023
%t With[{m = 20}, Range[0, m]! * CoefficientList[Series[1/(1 - x - 2*x^2 - 3*x^3), {x, 0, m}], x]] (* _Amiram Eldar_, Oct 30 2023 *)
%o (PARI) my(x='x+O('x^25)); Vec(serlaplace(1/(1-x-2*x^2-3*x^3))) \\ _Michel Marcus_, Oct 30 2023
%Y Cf. A000073, A000142, A052567, A052585, A101822, A189886, A364324.
%K nonn
%O 0,3
%A _Enrique Navarrete_, Oct 29 2023