Second order Eulerian numbers

A companion to A340556

Peter Luschny, Feb. 2021

You can download this SageMath notebook from GitHub.

Indexing the second-order Eulerian numbers comes in three flavors: A008517 (following Riordan and Comtet), A201637 (following Graham, Knuth, and Patashnik) and A340556, extending the definition of Gessel and Stanley. Note that we use the definition A340556.

Notations

$\left\langle\!\! \left\langle n\atop k\right\rangle\!\! \right\rangle $ A340556 second-order Eulerian numbers

$\left\langle \! \left\langle x \right\rangle \! \right\rangle _n $ second-order Eulerian polynomials

$\left[ n\atop k\right]$ A132393 Stirling cycle numbers (unsigned first kind)

$\left\{ n\atop k \right\} $ A048993 Stirling set numbers (signed second kind)

$\left[ \! \left[ n\atop k\right] \! \right] $ (A008306, A106828) associated Stirling cycle numbers (see def. below)

$\left\{ \!\! \left\{ n\atop k\right\} \!\! \right\} $ (A008299, A137375) associated Stirling set numbers (see def. below)

$\operatorname{W}_{n,k}$ A269939 Ward numbers

$\operatorname{W}_{n}(x)$ Ward polynomials

Formula 1

$$\left\langle\!\! \left\langle n\atop k\right\rangle\!\! \right\rangle = 0 \ \ (n \lt 0),\ \left\langle\!\! \left\langle n\atop k\right\rangle\!\! \right\rangle = k^n \ \ (k = 0), \\ \left\langle\!\! \left\langle n\atop k\right\rangle\!\! \right\rangle \ = \ k\, \left\langle\!\! \left\langle n-1\atop k\right\rangle\!\! \right\rangle + (2n-k)\, \left\langle\!\! \left\langle n-1\atop k-1\right\rangle\!\! \right\rangle $$

@cached_function
def E2(n, k):
    if n < 0: return 0
    if k == 0: return k^n
    return k * E2(n - 1, k) + (2*n - k) * E2(n - 1, k - 1)


for n in range(10):
    print([E2(n, k) for k in (0..n)])

[1]
[0, 1]
[0, 1, 2]
[0, 1, 8, 6]
[0, 1, 22, 58, 24]
[0, 1, 52, 328, 444, 120]
[0, 1, 114, 1452, 4400, 3708, 720]
[0, 1, 240, 5610, 32120, 58140, 33984, 5040]
[0, 1, 494, 19950, 195800, 644020, 785304, 341136, 40320]
[0, 1, 1004, 67260, 1062500, 5765500, 12440064, 11026296, 3733920, 362880]

Formula 2

$$\left\langle \! \left\langle x \right\rangle \! \right\rangle _n = \sum_{k=0}^{n} \left\langle\!\! \left\langle n\atop k\right\rangle\!\! \right\rangle x^k $$

# row polynomials:
def E2poly(n):
    return sum(E2(n, k) * x^k for k in (0..n)) 


for n in range(8):
    print(E2poly(n))

1
x
2*x^2 + x
6*x^3 + 8*x^2 + x
24*x^4 + 58*x^3 + 22*x^2 + x
120*x^5 + 444*x^4 + 328*x^3 + 52*x^2 + x
720*x^6 + 3708*x^5 + 4400*x^4 + 1452*x^3 + 114*x^2 + x
5040*x^7 + 33984*x^6 + 58140*x^5 + 32120*x^4 + 5610*x^3 + 240*x^2 + x

Formula 3

$$ \left\langle \! \left\langle x \right\rangle \! \right\rangle _n = x\,(x-1)^{2n} \, \left(\frac{ \left\langle \! \left\langle x \right\rangle \! \right\rangle _{n-1}}{(1-x)^{2n-1}} \right)' \quad (n \ge 1). $$

@cached_function
def e2poly(n):
    R.<x> = PolynomialRing(ZZ)
    if n == 0: return R(1)
    return R(expand(x * (x - 1)^(2 * n) * derivative((1 - x)^(1 - 2 * n) * e2poly(n - 1))))


for n in range(8):
    print(e2poly(n))
    
for n in (0..7):
    print(e2poly(n).coefficients(sparse=False))    

1
x
2*x^2 + x
6*x^3 + 8*x^2 + x
24*x^4 + 58*x^3 + 22*x^2 + x
120*x^5 + 444*x^4 + 328*x^3 + 52*x^2 + x
720*x^6 + 3708*x^5 + 4400*x^4 + 1452*x^3 + 114*x^2 + x
5040*x^7 + 33984*x^6 + 58140*x^5 + 32120*x^4 + 5610*x^3 + 240*x^2 + x
[1]
[0, 1]
[0, 1, 2]
[0, 1, 8, 6]
[0, 1, 22, 58, 24]
[0, 1, 52, 328, 444, 120]
[0, 1, 114, 1452, 4400, 3708, 720]
[0, 1, 240, 5610, 32120, 58140, 33984, 5040]

# plot([e2poly(n) for n in range(6)], color=['blue','red'], 
#     legend_label='E2poly', ymin=-1, ymax=2)

Formula 4

$$ \left\langle\! \left\langle x \right\rangle\! \right\rangle_n = (n+1)!\, (1-t)^{2n + 1}\, [x^{n+1}]\, \operatorname{Reversion}_{x}(x + t - t \exp(x)) \quad (n \ge 0)$$

# The compositional inverse with respect to x of x + t - t*exp(x).

# Mathematica:
# S := InverseSeries[Series[x + t - t Exp[x], {x, 0, 7}], x];
# Table[(n+1)! (1-t)^(2n+1) Coefficient[S, x, n+1], {n, 0, 6}] 

Formula 5

$$\operatorname{Egf}(n) = (1 - x)^{2n + 1} \sum_{k=0}^{n} \frac{(x\,\exp(-x))^{k}\ k^{n+k}}{k!} \\ \left\langle\!\! \left\langle n\atop k\right\rangle\!\! \right\rangle = [x^k]\ \operatorname{Egf}(n) \quad (n > 0) $$

var('x, k')

def Egf(n):
    if n == 0: return 1
    s = sum((x * exp(-x))^k * k^(n + k) / factorial(k) for k in (0..n))
    return (1 - x)^(2 * n + 1) * s


def E2row(n):
    return taylor(Egf(n), x, 0, n).list()

for n in (0..7):
    print(E2row(n))    

[1]
[0, 1]
[0, 1, 2]
[0, 1, 8, 6]
[0, 1, 22, 58, 24]
[0, 1, 52, 328, 444, 120]
[0, 1, 114, 1452, 4400, 3708, 720]
[0, 1, 240, 5610, 32120, 58140, 33984, 5040]

Formula 6

$$\left\langle\!\! \left\langle n\atop k\right\rangle\!\! \right\rangle \ = \ \sum \limits_{j\,=\,0}^{n-k} (-1)^j \binom{2n+1}{j}\genfrac[]{0pt}{}{n+m}{m} $$

where $m = n-k-j+1$.

# We denote alternative implementations of E2(n,k) by e2(n,k) (so we can test E2 = e2).

def e2(n, k):
    return sum((-1)^j * binomial(2 * n + 1, j) * stirling_number1(2 * n - k - j + 1, n - k - j + 1) 
           for j in (0..n-k))


for n in range(10):
    print([e2(n, k) for k in (0..n)])
    #print([E2(n, k) for k in (0..n)])

[1]
[0, 1]
[0, 1, 2]
[0, 1, 8, 6]
[0, 1, 22, 58, 24]
[0, 1, 52, 328, 444, 120]
[0, 1, 114, 1452, 4400, 3708, 720]
[0, 1, 240, 5610, 32120, 58140, 33984, 5040]
[0, 1, 494, 19950, 195800, 644020, 785304, 341136, 40320]
[0, 1, 1004, 67260, 1062500, 5765500, 12440064, 11026296, 3733920, 362880]

Formula 7

$$\left\langle\!\! \left\langle n\atop k\right\rangle\!\! \right\rangle \ = \ \sum \limits_{j\,=\,0}^{k} (-1)^{k-j}\binom{2n+1}{k-j} \genfrac\{\}{0pt}{}{n+j}{j} $$

def e2(n, k):
    return sum((-1)^(k - j) * binomial(2*n + 1, k - j)*stirling_number2(n + j, j)
           for j in (0..k))


for n in range(10):
    print([e2(n, k) for k in (0..n)])

[1]
[0, 1]
[0, 1, 2]
[0, 1, 8, 6]
[0, 1, 22, 58, 24]
[0, 1, 52, 328, 444, 120]
[0, 1, 114, 1452, 4400, 3708, 720]
[0, 1, 240, 5610, 32120, 58140, 33984, 5040]
[0, 1, 494, 19950, 195800, 644020, 785304, 341136, 40320]
[0, 1, 1004, 67260, 1062500, 5765500, 12440064, 11026296, 3733920, 362880]

Formula 8

$$\sum \limits_{j\,=\,0}^{k} (-1)^{k-j}\binom{2n+1}{k-j}{n+j \brace j} = \sum \limits_{j\,=\,0}^{n-k} (-1)^j \binom{2n+1}{j}{n+m \brack m}, $$

where $m = n-k-j+1$.

A proof of this identity was given by Marko Riedel at MSE.

Formula 9

$$ \left[ n\atop n-k\right] = \sum_{j=0}^k \binom{n + j - 1}{2k} \left\langle\!\! \left\langle k\atop j\right\rangle\!\! \right\rangle $$

def rhs(n, k): 
    return sum(E2(k, j) * binomial(n + j - 1, 2 * k) for j in (0..k))


for n in range(10):
    print([rhs(n, k) for k in (0..n)])
    #print([stirling_number1(n, n - k) for k in (0..n)])

[1]
[1, 0]
[1, 1, 0]
[1, 3, 2, 0]
[1, 6, 11, 6, 0]
[1, 10, 35, 50, 24, 0]
[1, 15, 85, 225, 274, 120, 0]
[1, 21, 175, 735, 1624, 1764, 720, 0]
[1, 28, 322, 1960, 6769, 13132, 13068, 5040, 0]
[1, 36, 546, 4536, 22449, 67284, 118124, 109584, 40320, 0]

Formula 10

$$ \genfrac\{\}{0pt}{}{n}{n-k} = \sum_{j=0}^k \binom{n + k - j}{2k} \left\langle\!\! \left\langle k\atop j\right\rangle\!\! \right\rangle $$

def rhs(n, k):
    return sum(binomial(n + k - j, 2*k) * E2(k, j) for j in (0..k))


for n in range(10):
    print([rhs(n, k) for k in (0..n)])
    #print([stirling_number2(n, n - k) for k in (0..n)])

[1]
[1, 0]
[1, 1, 0]
[1, 3, 1, 0]
[1, 6, 7, 1, 0]
[1, 10, 25, 15, 1, 0]
[1, 15, 65, 90, 31, 1, 0]
[1, 21, 140, 350, 301, 63, 1, 0]
[1, 28, 266, 1050, 1701, 966, 127, 1, 0]
[1, 36, 462, 2646, 6951, 7770, 3025, 255, 1, 0]

Formula 11

$$ \left| n\atop k\right| = \sum_{j=0}^k \left( \binom{n + j - 1}{2k} - \binom{n + k - j}{2k} \right) \left\langle\!\! \left\langle k\atop j\right\rangle\!\! \right\rangle $$

$$ \left[ n\atop k\right] - \left\{ n\atop k\right\} = \left| n\atop n- k\right| $$

A proof of this identity was given by René Gy and Marko Riedel at MSE. See also A341102.

Formula 12

$$ \operatorname{W}(n, 0) = 0^n, \operatorname{W}(n, k) = 0 \ \ (k < 0 \text{ or } k > n),$$

$$ \operatorname{W}(n, k) = k\,\operatorname{W}(n-1, k) + (n + k - 1)\,\operatorname{W}(n-1, k-1)$$

# The Ward numbers
@cached_function
def W(n, k):
    if k == 0 and n == 0: return 1
    if k <= 0 or  k >  n: return 0
    return k * W(n-1, k) + (n + k - 1) * W(n-1, k-1)


for n in (0..7):
    print([W(n, k) for k in (0..n)])

[1]
[0, 1]
[0, 1, 3]
[0, 1, 10, 15]
[0, 1, 25, 105, 105]
[0, 1, 56, 490, 1260, 945]
[0, 1, 119, 1918, 9450, 17325, 10395]
[0, 1, 246, 6825, 56980, 190575, 270270, 135135]

Formula 13

$$ \operatorname{W}_{n, k} = \sum_{j=0}^k (-1)^{j+k} \binom{n+k}{n+j} \left\{ n+j \atop j\right \} $$

def w(n, k):
    return sum((-1)^(j+k) * binomial(n + k, n + j) * stirling_number2(n + j, j)
           for j in (0..k))


for n in (0..7):
    print([w(n, k) for k in (0..n)])

[1]
[0, 1]
[0, 1, 3]
[0, 1, 10, 15]
[0, 1, 25, 105, 105]
[0, 1, 56, 490, 1260, 945]
[0, 1, 119, 1918, 9450, 17325, 10395]
[0, 1, 246, 6825, 56980, 190575, 270270, 135135]

Formula 14

$$ \operatorname{W}_{n,k} = \sum_{j=0}^{k} \binom{n-j}{n-k} \left\langle\!\! \left\langle n\atop j\right\rangle\!\! \right\rangle $$

def W(n, k):
    return sum(binomial(n - j, n - k) * E2(n, j) for j in (0..k))


for n in range(8):
    print([W(n, k) for k in (0..n)])

[1]
[0, 1]
[0, 1, 3]
[0, 1, 10, 15]
[0, 1, 25, 105, 105]
[0, 1, 56, 490, 1260, 945]
[0, 1, 119, 1918, 9450, 17325, 10395]
[0, 1, 246, 6825, 56980, 190575, 270270, 135135]

Formula 15

$$ \sum_{j=0}^{k} \binom{n-j}{n-k} \left\langle\!\! \left\langle n\atop j\right\rangle\!\! \right\rangle \,=\, \sum_{j=0}^k (-1)^{j+k} \binom{n+k}{n+j} \left\{ n+j \atop j\right \} \quad ( n \ge 0) $$

See Marko Riedel's proof at MSE.

Formula 16

$$ \operatorname{W}_{n}(x) = \sum_{k=0}^n \operatorname{W}(n, k)\, x^k $$

# The Ward polynomials
def Wpoly(n):
    return sum(W(n, k) * x^k for k in (0..n))


for n in range(8):
    print(Wpoly(n))
    
for n in (0..7):
    print(Wpoly(n).coefficients(sparse=False))        

1
x
3*x^2 + x
15*x^3 + 10*x^2 + x
105*x^4 + 105*x^3 + 25*x^2 + x
945*x^5 + 1260*x^4 + 490*x^3 + 56*x^2 + x
10395*x^6 + 17325*x^5 + 9450*x^4 + 1918*x^3 + 119*x^2 + x
135135*x^7 + 270270*x^6 + 190575*x^5 + 56980*x^4 + 6825*x^3 + 246*x^2 + x
[1]
[0, 1]
[0, 1, 3]
[0, 1, 10, 15]
[0, 1, 25, 105, 105]
[0, 1, 56, 490, 1260, 945]
[0, 1, 119, 1918, 9450, 17325, 10395]
[0, 1, 246, 6825, 56980, 190575, 270270, 135135]

# plot([Wpoly(n) for n in range(6)], color=['blue','red'], 
#     legend_label='Wpoly', ymin=-1, ymax=2)

Formula 17

$$ \operatorname{W}_{n}(x) = (1 + x)^n \left\langle \! \left\langle \frac{x}{1+x} \right\rangle \! \right\rangle _n $$

def wpoly(n):
    return expand(factor((1 + x)^n * E2poly(n).subs(x = x/(1 + x)))) 


for n in range(8):
    print(wpoly(n).coefficients(sparse=False))        

[1]
[0, 1]
[0, 1, 3]
[0, 1, 10, 15]
[0, 1, 25, 105, 105]
[0, 1, 56, 490, 1260, 945]
[0, 1, 119, 1918, 9450, 17325, 10395]
[0, 1, 246, 6825, 56980, 190575, 270270, 135135]

Formula 18

$$ \left\langle \! \left\langle x \right\rangle \! \right\rangle _n = (1 - x)^n \operatorname{W}_{n}\left(\frac{x}{1-x}\right) $$

# Conversely:
def e2poly(n):
    return expand(factor((1 - x)^n * Wpoly(n).subs(x = x / (1 - x)))) 


for n in range(8):
    print(e2poly(n))
    
for n in (0..7):
    print(e2poly(n).coefficients(sparse=False))

1
x
2*x^2 + x
6*x^3 + 8*x^2 + x
24*x^4 + 58*x^3 + 22*x^2 + x
120*x^5 + 444*x^4 + 328*x^3 + 52*x^2 + x
720*x^6 + 3708*x^5 + 4400*x^4 + 1452*x^3 + 114*x^2 + x
5040*x^7 + 33984*x^6 + 58140*x^5 + 32120*x^4 + 5610*x^3 + 240*x^2 + x
[1]
[0, 1]
[0, 1, 2]
[0, 1, 8, 6]
[0, 1, 22, 58, 24]
[0, 1, 52, 328, 444, 120]
[0, 1, 114, 1452, 4400, 3708, 720]
[0, 1, 240, 5610, 32120, 58140, 33984, 5040]

Formula 19

$$ \left\langle\!\! \left\langle n\atop k\right\rangle\!\! \right\rangle = \sum_{j=0}^{k}(-1)^{k-j} \binom{n-j}{n-k} \operatorname{W}_{n,j} $$

def e2(n, k):
    return sum((-1)^(k - j) * W(n, j) * binomial(n - j, k - j) for j in (0..k))


for n in range(8):
    print([e2(n, k) for k in (0..n)])

[1]
[0, 1]
[0, 1, 2]
[0, 1, 8, 6]
[0, 1, 22, 58, 24]
[0, 1, 52, 328, 444, 120]
[0, 1, 114, 1452, 4400, 3708, 720]
[0, 1, 240, 5610, 32120, 58140, 33984, 5040]

Formula 20

$$ \left[ \! \left[ n\atop k\right] \! \right] = (n-1)\, \left[ \! \left[ n-1\atop k\right] \! \right] + (n-1) \, \left[ \! \left[ n-2\atop k-1\right] \! \right]$$

# [[n, k]] are the associated Stirling numbers of the first kind.

@cached_function
def AS1(n, k):
    if n < 0: return 0
    if k == 0: return k^n
    return (n - 1) * AS1(n - 1, k) + (n - 1) * AS1(n - 2, k - 1)


for n in range(8):
    print([AS1(n,k) for k in (0..n)])

[1]
[0, 0]
[0, 1, 0]
[0, 2, 0, 0]
[0, 6, 3, 0, 0]
[0, 24, 20, 0, 0, 0]
[0, 120, 130, 15, 0, 0, 0]
[0, 720, 924, 210, 0, 0, 0, 0]

Formula 21

$$ \left[ \! \left[ {n \atop k} \right] \! \right] = \sum_{j=0}^{n-k} \binom{j}{n-2k} \left\langle \!\!\! \left\langle {n-k \atop j+1} \right\rangle \!\!\! \right\rangle $$

def aS1(n, k):
    if k == 0: return k^n
    return sum(binomial(j, n - 2 * k) * E2(n - k, j + 1) for j in (0..n-k))


for n in range(8):
    print([aS1(n, k) for k in (0..n)])
    #print([AS1(n, k) for k in (0..n)])

[1]
[0, 0]
[0, 1, 0]
[0, 2, 0, 0]
[0, 6, 3, 0, 0]
[0, 24, 20, 0, 0, 0]
[0, 120, 130, 15, 0, 0, 0]
[0, 720, 924, 210, 0, 0, 0, 0]

Formula 22

$$ \left\langle \!\!\left\langle {n \atop k} \right\rangle \!\!\right\rangle = \sum_{j=0}^{n-k+1} (-1)^{n-k-j+1} \binom{n-j}{k-1} \left[\!\left[ {n+j \atop j} \right]\!\right] \quad (n \ge 1)$$

def e2(n, k): 
    if n == k: return factorial(n)
    return sum((-1)^(n - k - j + 1)*binomial(n - j, k - 1) * AS1(n + j, j) for j in (0..n - k + 1))


for n in range(8):
    print([e2(n, k) for k in (0..n)])

[1]
[0, 1]
[0, 1, 2]
[0, 1, 8, 6]
[0, 1, 22, 58, 24]
[0, 1, 52, 328, 444, 120]
[0, 1, 114, 1452, 4400, 3708, 720]
[0, 1, 240, 5610, 32120, 58140, 33984, 5040]

Formula 23

$$ \left\{ \!\! \left\{ n\atop k\right\} \!\! \right\} = k\, \left\{ \!\! \left\{ n-1\atop k\right\} \!\! \right\} + (n+1) \, \left\{ \!\! \left\{ n-2\atop k-1\right\} \!\! \right\}$$

# {{n, k}} are the associated Stirling numbers of the second kind.

@cached_function
def AS2(n, k):
    if n < 0: return 0
    if k == 0: return k^n
    return k * AS2(n - 1, k) + (n - 1) * AS2(n - 2, k - 1)


for n in range(8):
    print([AS2(n, k) for k in (0..n)])

[1]
[0, 0]
[0, 1, 0]
[0, 1, 0, 0]
[0, 1, 3, 0, 0]
[0, 1, 10, 0, 0, 0]
[0, 1, 25, 15, 0, 0, 0]
[0, 1, 56, 105, 0, 0, 0, 0]

Formula 24

$$ \left\{ \!\! \left\{ n\atop k\right\} \!\! \right\} = \sum_{j=0}^{n-k}\binom{j}{n-2k} \left\langle\!\!\!\genfrac<>{0pt}{}{n-k}{n-k-j}\!\!\!\right\rangle $$

def aS2(n, k):
    return sum(binomial(j, n - 2 * k) * E2(n - k, n - k - j) for j in (0..n-k))


for n in range(8):
    print([aS2(n, k) for k in (0..n)])

[1]
[0, 0]
[0, 1, 0]
[0, 1, 0, 0]
[0, 1, 3, 0, 0]
[0, 1, 10, 0, 0, 0]
[0, 1, 25, 15, 0, 0, 0]
[0, 1, 56, 105, 0, 0, 0, 0]

Formula 25

$$ \left\langle\!\! \left\langle n\atop k\right\rangle\!\! \right\rangle = \sum_{j=0}^{k} (-1)^{k-j} \binom{n-j}{k-j} \left\{ \!\! \left\{ n+j\atop j\right\} \!\! \right\} $$

def e2(n, k):
    return sum((-1)^(k - j) * binomial(n - j, k - j) * AS2(n + j, j) for j in (0..k))


for n in range(8):
    print([e2(n, k) for k in (0..n)])

[1]
[0, 1]
[0, 1, 2]
[0, 1, 8, 6]
[0, 1, 22, 58, 24]
[0, 1, 52, 328, 444, 120]
[0, 1, 114, 1452, 4400, 3708, 720]
[0, 1, 240, 5610, 32120, 58140, 33984, 5040]

Formula 26

$$ \sum_{j=0}^{k} (-1)^{k-j} \binom{n-j}{k-j} \left\{ \!\! \left\{ n+j\atop j\right\} \!\! \right\} = \sum_{j=0}^{n-k+1} (-1)^{n-k-j+1} \binom{n-j}{k-1} \left[\! \left[ n+j\atop j\right] \! \right] $$

A proof of this identity was given by Marko Riedel at MSE.

Formula 27

Let $W(x) \,=\, $ LambertW$(x)$, then

$$ 2^n \left\langle\!\!\left\langle \frac{1}{2} \right\rangle\!\!\right\rangle_n = n! \, [x^n]\, \frac{1}{{2 + 2\operatorname{W}(- \exp((x-1)\,/\,2)\,/\,2)}} .$$

def lhs(n):
    return 2^n * E2poly(n).substitute(x=1/2) 


print([lhs(n) for n in range(11)])

# Maple:
# ser := series((1/2)/(1 + LambertW((-1/2)*exp((x-1)/2))), x, 14):
# seq(n!*coeff(ser, x, n), n=0..10);

[1, 1, 4, 26, 236, 2752, 39208, 660032, 12818912, 282137824, 6939897856]

Formula 28

$$ 2^n \left\langle\!\!\left\langle - \frac{1}{2} \right\rangle\!\!\right\rangle_n = (n + 1)!\, [x^{n+1}]\, \text{Reversion}\left(\frac{6x + \exp(3x) - 1}{9}\right) \quad (n \ge 0). $$

def lhs(n):
    return 2^n * E2poly(n).substitute(x=-1/2) 


def rhs(n):
    R.<x> = PowerSeriesRing(QQ, default_prec=2*n)
    f = R((6 * x + exp(3 * x) - 1) / 9)
    p = f.pade(n + 1, 0)  
    C = R(p).reverse(precision=n+1).shift(-1).list()
    return [c * factorial(k + 1) for k,c in enumerate(C)]


print([lhs(n) for n in range(11)])
print(rhs(11))

# Mathematica:
# f[x_] := (6 x + Exp[3 x] - 1)/9;
# S := InverseSeries[Series[f[x], {x, 0, R}], x];
# Table[(n + 1)!  Coefficient[S, x, n + 1], {n, 0, 10}]   

[1, -1, 0, 6, -12, -144, 1080, 5184, -127008, 95904, 19077120]
[1, -1, 0, 6, -12, -144, 1080, 5184, -127008, 95904, 19077120]

Formula 29

$$ \, 2^n \left\langle\!\!\left\langle - \frac{1}{2} \right\rangle\!\!\right\rangle_n \, = \, 3^n W_n\left(-\frac13\right) \quad ( n \ge 0) $$

def lhs(n):
    return 3^n * Wpoly(n).substitute(x=-1/3) 

print([lhs(n) for n in range(11)])

[1, -1, 0, 6, -12, -144, 1080, 5184, -127008, 95904, 19077120]

Formula 30

$$ m^{m+n} \, = \, m! \, [x^m] \frac{\left\langle \! \left\langle \,\operatorname{T}(x) \, \right\rangle \! \right\rangle _n} {(1-\operatorname{T}(x))^{2n+1}} \quad ( n \ge 0) $$

Here $\operatorname{T}(z) = - \operatorname{W}(-z)$ denotes Euler's tree function, and $\operatorname{W}$ is the principal branch of Lambert's function (A000169). See also Marko Riedel's proof on MSE.

def lhs(n, m):
    return m^(m + n)


for n in range(5):
    print([lhs(n, k) for k in (0..6)])

    
def rhs(n, L):
    R.<x> = PowerSeriesRing(QQ, L + (n == 0))
    T = R(sum(m^(m - 1) * x^m / factorial(m) for m in (1..L)))
    S = R(E2poly(n)).substitute(x = T) / (1 - T)^(2 * n + 1) 
    return [c * factorial(m) for m,c in enumerate(S.list())]


print()
for n in range(5): print(rhs(n, 6))

    
# Maple:
# T := x -> -LambertW(-x);
# ET := n -> E2poly(n, T(x)) / (1 - T(x))^(2*n+1);
# for n from 0 to 6 do
#    ser := series(ET(n), x, 16):
#    print(seq(k!*coeff(ser, x, k), k=0..9)) od:

[1, 1, 4, 27, 256, 3125, 46656]
[0, 1, 8, 81, 1024, 15625, 279936]
[0, 1, 16, 243, 4096, 78125, 1679616]
[0, 1, 32, 729, 16384, 390625, 10077696]
[0, 1, 64, 2187, 65536, 1953125, 60466176]

[1, 1, 4, 27, 256, 3125, 46656]
[0, 1, 8, 81, 1024, 15625, 279936]
[0, 1, 16, 243, 4096, 78125, 1679616]
[0, 1, 32, 729, 16384, 390625, 10077696]
[0, 1, 64, 2187, 65536, 1953125, 60466176]

Discussions on MSE

Topics tagged

Second-order Eulerian numbers, Lambert's W function, and Schröder's fourth problem

A Stirling number identity representing the second-order Eulerian numbers

A recurrence of the second-order Eulerian polynomials

The value of the second-order Eulerian polynomials at x = -1/2

An identity related to the second-order Eulerian numbers

Difference of the Stirling cycle numbers and the Stirling set numbers

An associated Stirling number identity related to the second-order Eulerian numbers