(PARI) a(n) = sum(k=0, n, binomial(n, k) * k^(n-k) * 3^(k*(n-k)) )
for(n=0, 25, print1(a(n), ", "))
(PARI) a(n)=n!*polcoeff(sum(k=0, n, exp(k*3^k*x +x*O(x^n))*x^k/k!), n)
for(n=0, 25, print1(a(n), ", "))
(PARI) a(n)=polcoeff(sum(k=0, n, x^k/(1-k*3^k*x +x*O(x^n))^(k+1)), n)
for(n=0, 25, print1(a(n), ", "))