This site is supported by donations to The OEIS Foundation.

User:Peter Luschny/EulerianPolynomialsGeneralized

From OeisWiki

Jump to: navigation, search

Generalized Eulerian Polynomials

The Eulerian polynomials were introduced by Leonhard Euler in his Remarques sur un beau rapport entre les séries des puissances tant directes que réciproques in 1749 (first printed in 1765) where he describes a method of computing values of the zeta function at negative integers by a precursor of Abel's theorem applied to a divergent series. The Eulerian polynomials should not be confused with the Euler polynomials.

In this second part we look at an analytic generalization of the polynomials. The first part of this post can be found here. The full post (rendered with MathJax) is on the author's homepage.

Contents

The connection with cardinal B-splines

The cardinal B-spline of the first order b1(n) is the characteristic function of the unit interval. The cardinal B-spline of order n > 1 is  b_n(x) = \int_{0}^{1} b_{n-1}(x-t) \, dt . Then for n > 0

 A_n(x) = n! \, \sum_{k=0}^{n-1} b_{n+1}(k+1)x^k .

This representation of the Eulerian polynomials suggests to look also at the midpoint Eulerian polynomials

 M_n(x) = 2^{n} n! \, \sum_{k=0}^{n} b_{n+1}(k+1/2)x^k .

The midpoint Eulerian polynomials

Generating function

The midpoint Eulerian polynomials are defined by the generating function

 \sum_{n=0}^{\infty} \frac{M_{n}(t)}{(t-1)^n} \, \frac{x^n}{2^n n!} = \frac{t-1}{t-e^x} e^{x/2},

which is the counterpart to the generating function of the standard Eulerian polynomials

 \sum_{n=0}^{\infty} \frac{A_{n}(t)}{(t-1)^n} \, \frac{x^n}{n!} = \frac{t-1}{t-e^x}.

Exponential generating function

 \sum_{n=0}^{\infty} M_{n}(x)\ \frac{t^n}{n!} = \frac{(1-x)\exp((1-x)t)}{1-x\exp(2(1-x)t)}.

Recurrence relation

The midpoint Eulerian polynomials can be computed by recurrence:

M0(t) = 1,
 M_{n}(t) = 2t(1-t)M'_{n-1}(t)+M_{n-1}(t)(1+(2n-1)t) \quad (n \ge 1)\,.

Expansion

The expansion analogous to Euler's given above is

 \sum_{j \ge 0} x^j(2j+1)^n = \frac{M_n(x)}{(1-x)^{n+1}}.

For instance we get for   n\,=\,0,1,2:

 1 + x + x^2 + x^3 + ... = \frac{1}{1-x},
 1 + 3x + 5x^2 + 7x^3 + ... = \frac{1+x}{(1-x)^2},
 1 + 3^2 x + 5^2 x^2 + 7^2 x^3 + ... = \frac{1+6x+x^2}{(1-x)^3}.

Worpitzky-type identity

 (2x+1)^{n} = \sum_{0 \le k \le n } \binom{x+k}{n} [x^k] M_n(x).

Roots of the polynomials

Mn(x) has n zeros which are simple and negative.

The midpoint Eulerian numbers

The coefficients of the midpoint Eulerian polynomials are the midpoint Eulerian numbers Mn,k

 M_{n}(t) = \sum_{k=0}^{n} M_{n,k}\ t^{k}.
A060187        A000165
Mn,k 01234row sum
010000 1
111000 2
216100 8
31232310 48
4176230761 384

Note that the offset in our setup is (n,k) = (0,0) which differs from the offset assumed in some of the formulas in A060187.

The combinatorial interpretation

Let Bn denote the set of signed permutations of

 I = \{ \pm i \ : \ 1 \le i \le n \}

such that p( − i) = − p(i) for all  i \in I . The descent number of p is defined as

 des(p) = \text{card} \{i \in [n]: p(i-1) > p(i)\}

where p(0) = 0. Then

 M_n(x) = \sum_{p \in B_n} x^{des(p)}.

The table below illustrates this representation for the case n = 2.

Signed permutations B2
and descents
pdespdes
-2, -1, 1, 2 01, -2, 2, -1 1
-2, 1, -1, 2 11, 2, -2, -1 1
-1, -2, 2, 1 12, -1, 1, -2 1
-1, 2, -2, 1 12, 1, -1, -2 2

The connection with the Lerch transcendent

Since the Lerch transcendent  \Phi(z,s,a) = \sum_{k=0}^{\infty} \frac{z^k}{(a+k)^s} (the analytic continuation is always implied) is related to the polylogarithm by   \operatorname{Li}_s(z) = z\Phi(z,s,a)  the formula given above shows

 \Phi(z,-n,1) = { {A_n(z)} \over (1-z)^{n+1}} \qquad (n \ge 0) ~.

Similarly one finds for the midpoint Eulerian polynomials

 2^n \Phi(z,-n,1/2) = { {M_n(z)} \over (1-z)^{n+1}} \qquad (n \ge 0) ~.

The Eulerian polynomials generalized

The formulas in the last section suggest to introduce the following generalization of the Eulerian polynomials An,k(z).

 k^n \Phi(z,-n,1/k) = { {A_{n,k}(z)} \over (1-z)^{n+1}} \qquad (n \ge 0, k \ge 1, z \ne 1) \qquad \text{(*)}

Thus An,1(z) = An(z) and An,2(z) = Mn(z).

A225117
 n  An,3(x)
0 1
1 2x + 1
2 4x2 + 13x + 1
3 8x3 + 93x2 + 60x + 1
4 16x4 + 545x3 + 1131x2 + 251x + 1
A225118
 n  An,4(x)
0 1
1 3x + 1
2 9x2 + 22x + 1
3 27x3 + 235x2 + 121x + 1
4 81x4 + 1996x3 + 3446x2 + 620x + 1

Plugging z = 1 into the left hand side of (*) times (1 − z)n + 1 would set An,k(1) = 0. On the other hand the row sums of the coefficients of An,k(z) are knn! (the k-factorial numbers, A000142, A000165, A032031,…). But taking the limit will reconcile the two values

 A_{n,k}(1) = \lim_{z \rightarrow 1} \ k^n (1-z)^{n+1} \Phi(z,-n,1/k) = k^n n! .

The most interesting point to evaluate these polynomials at is z = − 1 (this is what Euler did in the classical case).

An,1( − 1) A155585 1, 1, 0, -2, 0, 16, 0, -272, 0, 7936, 0,
An,2( − 1) A002436 1, 0, -4, 0, 80, 0, -3904, 0, 354560, 0,
An,3( − 1) A000810 1, -1, -8, 26, 352, -1936, -38528, 297296,
An,4( − 1) A079858 1, -2, -12, 88, 912, -11552, -176832,

(Note that in OEIS sequences are written zero-suppressed if every other term is 0 and that we refer to a sequence in the database even if it has a different sign pattern.)

If  0 \le n  is integer and 0 < a is real then [1]

 \Phi(z,-n,a) = a^n + \sum_{j=0}^{n} \binom{n}{j} \operatorname{Li}_{-j}(z) a^{n-j}.

Thus with k = 1 / a for some integer k > 0 the generalized Eulerian polynomials can be written

 A_{n,k}(z) = \left(1-z \right)^{n+1} \left(1 + \sum_{j=0}^{n} \binom{n}{j} \operatorname{Li}_{-j}(z) k^{j} \right) \qquad (z \ne 1).

Here the index k is used for integer values of k but it can also be seen as an arbitrary positive real. For example we get

2nAn,1 / 2( − 1) A119881 1, 3, 8, 18, 32, 48, 128, 528, 512,
3nAn,1 / 3( − 1) A225116 1, 5, 24, 110, 480, 2000, 8064, 32240,
2nAn,1 / 2(2) A052841 1, 0, 2, 6, 38, 270, 2342, 23646,

Assorted values of the polynomials

A(n, k, I) *
(1+I)^(1-n)
k = 1 k = 2 k = 3 k = 4
 Re A000111 A001586 A007286 A006873
 Im A000007 A001586 A007289 A225109
A(n, k, I) *
(1-I)^(1-n)
k = 1 k = 2 k = 3 k = 4
 Re A155585 A188458 A000810 A000813
A156205
 Im A122045 A212435 A225147 A156201

Among other numbers we see the Euler (or up/down) numbers, the generalized Euler numbers (Springer numbers), the number of alternating k-signed permutations and two recent additions triggered by this exploration. (As always the references are modulo sign and suppressed zeros.)

How to compute Eulerian polynomials.

The definitions given above are not computer friendly. Maple for example makes errors when computing the Lerch transcendent (as far as I know in all versions throughout the years 1998--2012) and gives wrong values if formula (*) is used. A computation based on the polylogarithm is much better, still the polylogarithm is not a trivial function to implement.

To the rescue come the cardinal splines (see the Schoenberg reference; Schoenberg elucidates the cases k = 1 and k = 2 in detail).

Cardinal B-splines are easy to compute and the Eulerian polynomials can be based on them. Despite their modern name B-splines were first investigated by Nikolai Lobachevsky in the nineteenth century.

@CachedFunction    
def B(n, x):   # Cardinal B-splines
    if n == 1: return 0 if (x < 0) or (x >= 1) else 1 
    return (x/(n-1))*B(n-1, x)+((n-x)/(n-1))*B(n-1, x-1)

def EulerianPolynomial(n, k, x):
    if n == 0: return 1
    return k^n*factorial(n)*add(B(n+1, m+1/k)*x^m for m in (0..n))

For example:

[EulerianPolynomial(n, 4, x) for n in (0..5)] 
[imag(-(1-I)^(1-n)*EulerianPolynomial(n,4,I)) for n in (0..17)] # A156201

No library functions besides the power and factorial function are used. And the factor k^n \, n! can be integrated into the recurrence.

@CachedFunction    
def EB(n, k, x):
    if n == 1: return 0 if (x < 0) or (x >= 1) else 1 
    return k*x*EB(n-1, k, x) + k*(n-x)*EB(n-1, k, x-1)
    
def EulerianPolynomial(n, k, x):
    if n == 0: return 1
    return add(EB(n+1, k, m+1/k)*x^m for m in (0..n))

Also the rational parameter m + 1 / k can be eliminated.

@CachedFunction    
def BB(n, k, x):  
    if n == 1: return 0 if (x < 0) or (x >= k) else 1 
    return x*BB(n-1, k, x) + (n*k-x)*BB(n-1, k, x-k) 

def EulerianPolynomial(n, k, x):
    if n == 0: return 1
    return add(BB(n+1, k, k*m+1)*x^m for m in (0..n))

Based on this recurrence a closed form can also be given if we use the signum function defined  \operatorname{sgn}(x-y) = -1 \text{ if } x < y; 0 \text{ if } x = y; 1 \text{ if } x > y. Then

def EulerianPolynomial(n, k, x) :
    if n == 0: return 1
    S = lambda m: add((-1)^j*binomial(n+1,j)*(k*(m-j)+1)^n*sgn(k*(m-j)+1) 
        for j in (0..n+1))
    return add(S(m)*x^m for m in (0..n))/2

Finally calling the coefficients of the (generalized) Eulerian polynomials (generalized) Eulerian numbers we arrive at the recurrence:

@CachedFunction
def EulerianNumber(n, k, m) :
    if n == 0: return 1 if m == 0 else 0
    return (k*(n-m+1)-1)*EulerianNumber(n-1,k,m-1)
               + (m*k+1)*EulerianNumber(n-1,k,m)

def EulerianPolynomial(n, k, x):
    return add(EulerianNumber(n, k, m)*x^m for m in (0..n))    

The generalized Eulerian polynomials can also be computed by a direct recurrence:

A0,k(t) = 1,
 A_{n,k}(t) = kt(1-t)A'_{n-1,k}(t)+A_{n-1,k}(t)(1+(kn-1)t) \quad (n \ge 1)\,.

Generating function

The generating function for the generalized Eulerian polynomials An,k(x) is

 \operatorname{gf}(n; k) = k^n n! \,  \left( 1 /x - 1 \right) ^{n+1} \, [t^n] \, x e^{\frac {tx}{k} } \left( 1-x e^{tx} \right) ^{-1}.

Here [tn]f(t,x) is the coefficient of tn in f(t,x).

This generating function immediately suggests a more symmetric generalization with two parameters:

 \operatorname{gf}(n; k, m) = (km)^n n! \, \left( 1 /x - 1 \right) ^{n+1} \, [t^n] \, x e^{\frac {tx}{k} } \left( 1-x e^{\frac {tx}{m}} \right) ^{-1}.

Short table An,k(x)

\scriptstyle A_{0,k}(x) = 1
\scriptstyle A_{1,k}(x) = (k - 1)\,x + 1
\scriptstyle A_{2,k}(x) = (k^{2} - 2\,k + 1)\,x^{2} + (k^{2} + 2\,k - 2)\,x + 1
\scriptstyle A_{3,k}(x) = (k^{3} - 3\,k^{2} + 3\,k - 1)\,x^{3} + (4\,k^{3} - 6\,k + 3)\,x^{2} + (k^{3} + 3\,k^{2} + 3\,k - 3)\,x + 1
\scriptstyle A_{4,k}(x) = (k^{4} - 4\,k^{3} + 6\,k^{2} - 4\,k + 1)\,x^{4} + (11\,k^{4} - 12\,k^{3} - 6\,k^{2} + 12\,k - 4)\,x^{3} + \ldots
\scriptstyle \qquad\qquad \ldots \ (11\,k^{4} + 12\,k^{3} - 6\,k^{2} - 12\,k + 6)\,x^{2} + (k^{4} + 4\,k^{3} + 6\,k^{2} + 4\,k - 4)\,x + 1

Program

(Maple)
a := proc(n, m) local k; # Eulerian numbers
     add((-1)^k*binomial(n+1, k)*(m+1-k)^n, k=0..m) 
     end:
A := proc(n, x) local k; # Eulerian polynomials
     add(a(n, k)*x^k, k=0..n) 
     end:
ma := proc (n, m) local k; # Midpoint Eulerian numbers
     add((-1)^(m-k)*binomial(n+1, m-k)*(2*k+1)^n, k=0..m)
     end:
mr := proc(n, k) option remember; # Recursive mid. Eul.num.
     if n = 0 then if k=0 then 1 else 0 fi else
     (2*(n-k)+1)*mr(n-1, k-1) + (2*k+1)*mr(n-1, k) fi 
     end:
MA := proc(n, x) local k; # Midpoint Eulerian polynomials
     add(mr(n, k)*x^k, k=0..n)
     end:
B := proc(n, u) option remember; # Cardinal B-splines
     if n = 1 then if (u < 0) or (u >= 1) then 0 else 1 fi
     else (u/(n-1))*B(n-1, u)+((n-u)/(n-1))*B(n-1, u-1) fi
     end:
     
# Generalized Eulerian polynomials based on polylog.
EulerianPolynomial := proc(n, k, x) local j;
     if x = 1 then k^n*n! else 
     (1-x)^(1+n)*(1+add(binomial(n,j)*polylog(-j,x)*k^j,j=0..n)); 
     simplify(expand(%)) fi 
     end:

# Generalized Eulerian polynomials based on direct recurrence.     
EulerianPolynomials := proc(n, k, t)
     if n = 0 then 1 else
     k*t*(1-t)*diff(EulerianPolynomials(n-1,k,t),t)
     +EulerianPolynomials(n-1,k,t)*(1+(k*n-1)*t) fi; 
     sort(simplify(expand(%))) end:
   
# Generating function, general case
gf := proc(n, k) local f; 
     f := (x,t) -> x*exp(t*x/k)/(1-x*exp(t*x));
     series(f(x,t), t, n+2); ((1-x)/x)^(n+1)*k^n*n!*coeff(%, t, n); 
     collect(simplify(%), x) end:
     seq(print(gf(n, k)), n=0..4); # generates the table above
    
(Sage)
def a(n, m) :  # Eulerian numbers
    return add((-1)^k*binomial(n+1, k)*(m+1-k)^n for k in (0..m))
    
def A(n, x) :  # Eulerian polynomials
    return add(a(n, k)*x^k for k in (0..n))
    
def ma(n, m):  # Midpoint Eulerian numbers
    return add((-1)^(m-k)*binomial(n+1, m-k)*(2*k+1)^n for k in (0..m))

@CachedFunction
def mr(n, k) : # Recursive midpoint Eulerian numbers
    if n == 0: return 1 if k == 0 else 0
    return (2*(n-k)+1)*mr(n-1, k-1) + (2*k+1)*mr(n-1, k) 
        
def MA(n, x):  # Midpoint Eulerian polynomials
    return add(mr(n, k)*x^k for k in (0..n))

@CachedFunction    
def B(n, x):   # Cardinal B-splines
    if n == 1: return 0 if (x < 0) or (x >= 1) else 1 
    return (x/(n-1))*B(n-1, x)+((n-x)/(n-1))*B(n-1, x-1)
    
@CachedFunction    
def BB(n, k, x):  # Modified Cardinal B-splines
    if n == 1: return 0 if (x < 0) or (x >= k) else 1 
    return x*BB(n-1, k, x) + (n*k-x)*BB(n-1, k, x-k) 

# Generalized Eulerian polynomials based on recurrence.      
def EulerianPolynomial(n, k, x):
    if n == 0: return 1
    return add(BB(n+1, k, k*m+1)*x^m for m in (0..n))

# Generalized Eulerian polynomials based on direct recurrence.          
def EulerianPolynomials(n, k, t):
    if n == 0 : return 1
    return (k*t*(1-t)*diff(EulerianPolynomials(n-1,k,t),t)
     + EulerianPolynomials(n-1,k,t)*(1+(k*n-1)*t)).expand() 

# Generalized Eulerian polynomials based on polylog.    
from mpmath import *
mp.dps = 32; mp.pretty = True

def EulerianPolynomialLi(n, k, x):
    if x == 1: return k^n*factorial(n)
    return (1-x)^(1+n)*(1+add(binomial(n,j)*polylog(-j,x)*k^j 
           for j in (0..n)))       

Notes

  1. See the formula at Wolfram research.

References

Personal tools