%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