%I #10 May 11 2021 15:32:06

%S 0,0,1,8,64,540,4920,48720,524160,6108480,76809600,1037836800,

%T 15008716800,231437606400,3792255667200,65819609856000,

%U 1206547550208000,23297526540288000,472708591939584000,10055994967130112000,223826984752250880000,5202760944485744640000,126075414965721661440000,3179798058882852126720000,83346901966165164687360000,2267221868000212451328000000

%N a(n) = Sum_{k = 0..n} E1(n, k)*k^2, where E1 are the Eulerian numbers A173018.

%C The Eulerian transform of the squares.

%F a(n) = n! * [x^n] x^2*(-x^2 + x - 3)/(6*(x - 1)^3).

%F a(n) = Sum_{k=0..n} Sum_{j=0..k} (-1)^j*binomial(n + 1, j)*k^2*(k + 1 - j)^n.

%F a(n) = ((n - 3)*(n - 1)*(23*n - 44)*a(n-2) + ((159 - 7*n)*n - 286)*a(n-1))/(16*(n - 2)) for n >= 3.

%p a := n -> add(combinat[eulerian1](n, k)*k^2, k = 0..n):

%p # Recurrence:

%p a := proc(n) option remember; if n < 2 then 0 elif n = 2 then 1 else

%p ((n-3)*(n-1)*(23*n-44)*a(n-2) + ((159 - 7*n)*n - 286)*a(n-1))/(16*(n - 2)) fi end:

%p seq(a(n), n = 0..29);

%t a[n_] := Sum[Sum[(-1)^j Binomial[n + 1, j] k^2 (k + 1 - j)^n, {j,0,k}], {k,0,n}]; a[0] := 0; Table[a[n], {n, 0, 25}]

%o (SageMath)

%o def aList(len):

%o R.<x> = PowerSeriesRing(QQ, default_prec=len+2)

%o f = x^2*(-x^2 + x - 3)/(6*(x - 1)^3)

%o return f.egf_to_ogf().list()[:len]

%o print(aList(20))

%Y Transforms of the squares: A151881 (StirlingCycle), A033452 (StirlingSet), A105219 (Laguerre), A103194 (Lah), A065096 (SchröderBig), A083411 (Fubini), A141222 (Narayana), A000330 (Units A000012).

%Y Cf. A173018.

%K nonn

%O 0,4

%A _Peter Luschny_, May 11 2021