User:Peter Luschny/P-Transform

The P-Transform

♦ P-polynomials

We introduce a sequence-to-triangle transformation which we will call the partition-transformation (or short P-transformation) similar to the Bell transformation discussed in the previous blog post which however is designed for a different type of formal series. In principle we could completely parallelize the exposition of the Bell transformation and in some aspects we will do so.

We start, as any mathematical text, with definitions. Unlike mathematical texts we use the Sage language instead of traditional mathematical notation.

def PartialPtrans(n, f, norm = None):
    if n == 0: return [1]
    R = [0]*(n+1)
    for q in Partitions(n):
        p = q + [0]
        R[q[0]] += (-1)^q[0]*mul(binomial(p[j], p[j+1])*f(j+1)^p[j] 
                                 for j in range(len(q)))
    if norm == None: return R
    return [norm(n, k)*R[k] for k in (0..n)]
    
Ptrans = lambda n, f, norm = None: sum(PartialPtrans(n, f, norm))    

Here we introduce two notions, the partial P-transformation 'PartialPtrans' and the P-transformation 'Ptrans'. The name obviously derives from the main loop which runs over all partitions of n. Next we explain their meaning by giving examples.

The most general case is given by 'Ptrans(n, f)', n≥0 an integer and f the function which just supplies the symbolic variables x_n. It represents what we call the partition polynomials of degree n which depend (in complete analogy to the case of Bell polynomials) on an unbounded number of variables x₁, x₂, x₃, ....

f = lambda n: var('x'+str(n))
for n in range(5): print Ptrans(n, f)

  1, 
-x1, 
 x1^2 - x1*x2, 
-x1^3 + 2*x1^2*x2 - x1*x2*x3, 
 x1^4 - 3*x1^3*x2 + x1^2*x2^2 + 2*x1^2*x2*x3 - x1*x2*x3*x4. 

The rows of this triangle might also be called the 'complete P-polynomials' to distinguish them from the partial P-polynomials, which are listed in this triangle for 0 ≤ k ≤ n:

for n in range(5): print PartialPtrans(n, f)

[1]
[0, -x1]
[0, -x1*x2,       x1^2]
[0, -x1*x2*x3,    2*x1^2*x2,                -x1^3]
[0, -x1*x2*x3*x4, x1^2*x2^2 + 2*x1^2*x2*x3, -3*x1^3*x2, x1^4]   

Obviously the complete P-polynomial is just the sum of the corresponding partial P-polynomials in row n.

The next step is to call Ptrans(n, f) and PartialPtrans(n, f) with a concrete sequence f (defined as a function). It substitutes the variables x₁, x₂, x₃, ... by the values f(1), f(2), f(3), ... (in other words, it evaluates the (multivariate) P-polynomials at the values f(i)). The resulting (univariate) polynomial can in turn by evaluted at any value 'val'. We implement this as:

def PEvaluateAt(n, f, val = None, norm = None):
    P = PartialPtrans(n, f, norm)
    p = sum([x^k*g for k,g in enumerate(P)])
    if val == None: return p
    return p.substitute(x=val)

For linguistic convenience we add a macro for the corner case where no substitution takes place (val = None).

AsPolynomials = lambda n, f, norm=None: PEvaluateAt(n, f, norm = norm)

Now recall that 'lambda n: 1' denotes in the Sage language the constant sequence a(n) = 1 for all n and 'lambda n: n' the sequence a(n) = n written as an anonymous function. (Maple users would write this as 'n → n'.) In these cases we find:

for n in range(6): print AsPolynomials(n, lambda n: 1)
1, 
-x, 
x^2 - x, 
-x^3 + 2*x^2 - x, 
x^4 - 3*x^3 + 3*x^2 - x, 
-x^5 + 4*x^4 - 6*x^3 + 4*x^2 - x.

for n in range(6): print AsPolynomials(n, lambda n: n)    
1
-x
x^2 - 2*x
-x^3 + 4*x^2 - 6*x
x^4 - 6*x^3 + 16*x^2 - 24*x
-x^5 + 8*x^4 - 30*x^3 + 72*x^2 - 120*x

What remains to make the OEIS happy is now routine:

• extracting the coefficients of the polynomials,
• evaluating the polynomials at x=1,
• evaluating the polynomials at x=-1,
• evaluating the polynomials at x=1/2 and multiplying the result with 2^n,
• evaluating the polynomials at x=-1/2 and multiplying the result with 2^n.

This translates in Sage to:

print "As polynomials:"
for n in range(Len-2): print AsPolynomials(n, f, norm)
print "Coefficients:"    
for n in range(Len-2): print PartialPtransRec(n, f, norm)    
    
print "x=1:   ", [PEvaluateAt(n,f,   1,norm)     for n in range(Len)]
print "x=-1:  ", [PEvaluateAt(n,f,  -1,norm)     for n in range(Len)]
print "x=1/2: ", [PEvaluateAt(n,f, 1/2,norm)*2^n for n in range(Len)]
print "x=-1/2:", [PEvaluateAt(n,f,-1/2,norm)*2^n for n in range(Len)]

Let's look at our basic examples:

f = lambda n: 1
norm = None

As polynomials:
1
-x
x^2 - x
-x^3 + 2*x^2 - x
x^4 - 3*x^3 + 3*x^2 - x
-x^5 + 4*x^4 - 6*x^3 + 4*x^2 - x

Coefficients:
[[1], 
[0, -1], 
[0, -1, 1], 
[0, -1, 2, -1], 
[0, -1, 3, -3, 1], 
[0, -1, 4, -6, 4, -1]]  # A097805, A071919, 

x=1:    [1, -1, 0, 0, 0, 0, 0, 0, 0, 0]             # A154955 
x=-1:   [1, 1, 2, 4, 8, 16, 32, 64, 128, 256]       # A011782 
x=1/2:  [1, -1, -1, -1, -1, -1, -1, -1, -1, -1]     # A153881
x=-1/2: [1, 1, 3, 9, 27, 81, 243, 729, 2187, 6561]  # A133494

f = lambda n: n

As polynomials:
1
-x
x^2 - 2*x
-x^3 + 4*x^2 - 6*x
x^4 - 6*x^3 + 16*x^2 - 24*x
-x^5 + 8*x^4 - 30*x^3 + 72*x^2 - 120*x

Coefficients:
[[1], 
[0, -1], 
[0, -2, 1], 
[0, -6, 4, -1], 
[0, -24, 16, -6, 1], 
[0, -120, 72, -30, 8, -1]]  # A090238, A059369  

x=1:    [1, -1, -1, -3, -13, -71, -461, -3447, -29093, -273343]  # A167894, A003319 
x=-1:   [1, 1, 3, 11, 47, 231, 1303, 8431, 62391, 524495]        # A051296
x=1/2:  [1, -1, -3, -17, -139, -1449, -18131, -263233]
x=-1/2: [1, 1, 5, 33, 269, 2633, 30421, 408945, 6307549]

We have now discussed almost all the basic definitions concerning the P-transformation. But the reader certainly noticed a third parameter in the list of parameters of the functions PartialPtrans and Ptrans, the optional parameter 'norm' which we have not used until now. 'norm' stands for normalization and can be thought as an additional factor applied to the partial P-transformations the purpose of which is to assure that the result are integers. The next two examples demonstrate this.

♦ Representations of the Euler and Bernoulli numbers

In an attempt to convince the reader that the P-transformation is interesting we show how the Euler and the Bernoulli numbers are represented as P-transforms. These examples also show that when we are speaking of a sequence-to-triangle transformation we specifically allow rational sequences as input for the transformation.

Euler numbers

We define as the generator sequence for the Euler numbers the rational sequence a(n) = 1/((2n-1)(2n)) and normalize the result by multiplying with (2n)!. Then we evaluate.

f = lambda n: 1/((2*n-1)*(2*n))
norm = lambda n,k: factorial(2*n)
Len = 8

print "As polynomials:"
for n in range(Len-2): print AsPolynomials(n, f, norm)
print "Coefficients:"    
for n in range(Len-2): print PartialPtransRec(n, f, norm)    
    
print "x=1:   ", [PEvaluateAt(n,f,   1,norm)     for n in range(Len)]
print "x=-1:  ", [PEvaluateAt(n,f,  -1,norm)     for n in range(Len)]
print "x=1/2: ", [PEvaluateAt(n,f, 1/2,norm)*2^n for n in range(Len)]
print "x=-1/2:", [PEvaluateAt(n,f,-1/2,norm)*2^n for n in range(Len)]

And we get the results:

# As polynomials:
1
-x
6*x^2 - x
-90*x^3 + 30*x^2 - x
2520*x^4 - 1260*x^3 + 126*x^2 - x
-113400*x^5 + 75600*x^4 - 13230*x^3 + 510*x^2 - x

# Coefficients: A241171 Joffe's central differences of zero, signed.
[1]
[0, -1]
[0, -1, 6]
[0, -1, 30, -90]
[0, -1, 126, -1260, 2520]
[0, -1, 510, -13230, 75600, -113400]

# x=1: A000364 Euler numbers signed, A028296 Gudermannian
[1, -1, 5, -61, 1385, -50521, 2702765, -199360981]

# x=-1: A094088 E.g.f. 1/(2-cosh(x)) (even coefficients).
[1, 1, 7, 121, 3907, 202741, 15430207, 1619195761]

# x=1/2: A002105 Reduced tangent numbers, signed.
[1, -1, 4, -34, 496, -11056, 349504, -14873104]

# x=-1/2: missing in the database
[1, 1, 8, 154, 5552, 321616, 27325088, 3200979664]

A big advantage of using the P-transformation is that it makes immediately visible a whole set of relations and sequences at the same time with a minimum of technical effort. Even more so when we include the inverse relations which we will discuss in the next section.

For those who have not yet become accustomed to the Sage language here the representation of the Euler numbers that we have computed in mathematical notation:

$\textstyle E_{2n}=(2n)!\sum _{p\in P_{n}}(-1)^{p_{1}}{\binom {p_{1}}{p_{2}}}{\binom {p_{2}}{p_{3}}}\ldots {\binom {p_{n-1}}{p_{n}}}\left({\frac {1}{1\cdot 2}}\right)^{p_{1}}\left({\frac {1}{3\cdot 4}}\right)^{p_{2}}\ldots \left({\frac {1}{(2n-1)2n}}\right)^{p_{n}}$

Bernoulli numbers

Now let's consider the Bernoulli numbers. This case is a bit more complex as not only the input sequence is a rational sequence but also the associated P-polynomials have rational coefficients.

We define as the generator sequence for the Bernoulli numbers the rational sequence a(n) = 1/((2n)(2n+1)) and normalize the result by multiplying with (2n)!/(2-2^(2n)). Then we evaluate:

f = lambda n: 1/((2*n)*(2*n+1))
norm = lambda n,k: factorial(2*n)/(2-2^(2*n))

print "\n # As polynomials:"
for n in range(Len-2): print AsPolynomials(n, f, norm)
print "\n # Coefficients"    
for n in range(Len-2): print PartialPtransRec(n, f, norm)    
print "\n # x=1:  ", [PEvaluateAt(n,f, 1,norm) for n in range(Len)]
print "\n # x=-1: ", [PEvaluateAt(n,f,-1,norm) for n in range(Len)]

Resulting in:

# As polynomials:
1
1/6*x
-1/21*x^2 + 1/70*x
5/93*x^3 - 1/31*x^2 + 1/434*x
-140/1143*x^4 + 14/127*x^3 - 41/1905*x^2 + 1/2286*x
100/219*x^5 - 40/73*x^4 + 93/511*x^3 - 23/1533*x^2 + 1/11242*x

# Coefficients
[1]
[0, 1/6]
[0, 1/70, -1/21]
[0, 1/434, -1/31, 5/93]
[0, 1/2286, -41/1905, 14/127, -140/1143]
[0, 1/11242, -23/1533, 93/511, -40/73, 100/219]

# x=1: A000367 / A002445
[1, 1/6, -1/30, 1/42, -1/30, 5/66, -691/2730, 7/6]
# [1, 1, -1, 1, -1, 5, -691, 7]
# [1, 6, 30, 42, 30, 66, 2730, 6]

# x=-1:
[1, -1/6, -13/210, -115/1302, -2911/11430, -13509/11242]
# [1, -1, -13, -115, -2911, -13509, -46675903]
# [1, 6, 210, 1302, 11430, 11242, 5588310]

It is natural to look at the polynomials (or their coefficients) as a refinement or generalization of the value of these polynomials at x = 1. In this sense we found here a generalization of the Bernoulli numbers which is not yet in the OEIS. In fact we have found more: a family of polynomials with vanishing constant terms (for n>0) which sum at x = 1 to the Bernoulli numbers. This is quite different from the classical Bernoulli polynomials which have non-vanishing constant terms for even n.

The fine-structure of the Bernoulli numbers which we used in this calculation written in mathematical notation:

$\textstyle B_{2n}={\frac {(2n)!}{(2-2^{2n})}}\sum _{p\in P_{n}}(-1)^{p_{1}}{\binom {p_{1}}{p_{2}}}{\binom {p_{2}}{p_{3}}}\ldots {\binom {p_{n-1}}{p_{n}}}\left({\frac {1}{2\cdot 3}}\right)^{p_{1}}\left({\frac {1}{4\cdot 5}}\right)^{p_{2}}\ldots \left({\frac {1}{2n(2n+1)}}\right)^{p_{n}}$

♦ The inverse of the P-transform

The inverse of the P-transform is defined as the inverse of the lower triangular matrix of the partial P-polynomials.

def InversePartialPtrans(dim, f, norm = None):
    A = [PartialPtrans(n, f) for n in range(dim)]
    M = [[0 for k in (0..n)] for n in range(dim)]
    M[0][0] = 1
    for n in (1..dim-1):
        M[n][n] = 1/A[n][n] 
        for k in range(n-1, 0, -1):
            M[n][k] = expand(-sum(A[i][k]*M[n][i] 
                      for i in range(n, k, -1))/A[k][k])

    if norm == None: return M
    for n in (1..dim-1): 
        for k in (1..n): M[n][k] *= norm(n,k)
    return M

In its full generality the matrix of the inverse partial P-polynomials starts:

f = lambda n: var('x'+str(n))
for p in InversePartialPtrans(5, f): print p

[1],
[0, -1/x1],
[0, -x2/x1^2,                                  1/x1^2],
[0, -2*x2^2/x1^3 + x2*x3/x1^3,                 2*x2/x1^3,                -1/x1^3],
[0, -5*x2^3/x1^4+5*x2^2*x3/x1^4-x2*x3*x4/x1^4, 5*x2^2/x1^4-2*x2*x3/x1^4, -3*x2/x1^4, 1/x1^4]

We see that x₁ != 0 is a necessary condition for the existence of the inverse P-transform. Column 1 of this matrix will be called the P-inverse of the sequence x₁, x₂, x₃,... . The inverse P-matrix can be cleaned up a bit if we assume x₁ = 1 (which will often be the case in practice).

for P in InversePartialPtrans(5, f): print [p.substitute(x1=1) for p in P]

[1]
[0, -1]
[0, -x2,                            1]
[0, -2*x2^2 + x2*x3,                2*x2,             -1]
[0, -5*x2^3 + 5*x2^2*x3 - x2*x3*x4, 5*x2^2 - 2*x2*x3, -3*x2, 1]

For example for the identity on N we get:

M = InversePartialPtrans(8, lambda n: n) 
for n in range(8): print [m for m in M[n]]

[1],
[0,  -1],
[0,  -2,  1],
[0,  -2,  4,  -1],
[0,  -4,  8,  -6,  1],
[0,   4, 16, -18,  8,  -1],
[0, -48, 12, -44, 32, -10,  1],
[0, 336, 96, -72, 96, -50, 12, -1].

The triangle is row reversed A059370 with different signs. And -1, -2, -2, -4, 4, -48, 336, ... is the P-inverse sequence of the sequence a(n) = n which is A059372 with a different sign pattern and called "revert transform of factorials".

def InversePsequence(n, f, norm = None):
    M = InversePartialPtrans(n, f, norm)
    return [M[j][1] for j in (1..n-1)]
    
print InversePsequence(8, lambda n: n)
[-1, -2, -2, -4, 4, -48, 336]    

♦ Recurrence for the case x = 1

Evaluating the P-polynomials at x = 1 is the same as computing the sum of the coefficients, or the sum of the rows if the coefficients are given by a lower triangular matrix. This is what the function 'Pcoeffsum(Len, f)' accomplishes. The integer 'Len' requests a list of length Len and 'f' is the sequence to be transformed, given as a function. In other words, Pcoeffsum(Len, f) is equivalent to [PEvaluateAt(k, lambda n: n, 1) for k in range(Len)] but it uses a more efficient algorithm.

def Pcoeffsum(len, f, norm = None):
    R, C = [1], [1]+[0]*(len-1)
    for n in (1..len-1):
        for k in range(n, 0, -1):
            C[k] = C[k-1] * f(k)
        C[0] = -sum(C[k] for k in (1..n))
        if norm == None: R.append(C[0])
        else: R.append(C[0]*norm(n))     
    return R

An example call:

print Pcoeffsum(9, lambda n: n)
[1, -1, -1, -3, -13, -71, -461, -3447, -29093]
print [PEvaluateAt(n, lambda n: n, 1) for n in range(9)]
[1, -1, -1, -3, -13, -71, -461, -3447, -29093]

There are many sequences in the OEIS which can be computed efficiently with this algorithm. A small selection of these sequences is given below. Of course in practice the given f will be built directly into the transformation and some normalization has to be added. For example the Sage code for the cosecant numbers looks like this:

def A001896_list(len):
    R, C = [1], [1]+[0]*(len-1)
    for n in (1..len-1):
        for k in range(n, 0, -1):
            C[k] = C[k-1] / (8*k*(2*k+1))
        C[0] = -sum(C[k] for k in (1..n))
        R.append((C[0]*factorial(2*n)).numerator())
    return R

print A001896_list(8) 

Note that the MacLaurin numbers which are the coefficients appearing in MacLaurin's second summation formula can be computed similarly: one only has to change the sign of f, f(k) = -1/(8*k*(2*k+1)) (compare A036280).

Another example: Choosing f(k) = (k+1) we get the convolutory inverse of the factorial sequence A077607. But choosing f(k) = 1/(k+1) and normalizing with n! we get ... the Bernoulli numbers again:

def bernoulli_list(len):
    f, R, C = 1, [1], [1]+[0]*(len-1)
    for n in (1..len-1):
        f *= n
        for k in range(n, 0, -1):
            C[k] = C[k-1] / (k+1)
        C[0] = -sum(C[k] for k in (1..n))
        R.append(C[0]*f)
    return R

print bernoulli_list(17) 

The associated polynomials are displayed in the plot below.

P = plot([1/2*x^2 - 1/3*x, -3/4*x^3 + x^2 - 1/4*x, 3/2*x^4 - 3*x^3 + 5/3*x^2 
-1/5*x, -15/4*x^5 + 10*x^4 - 35/4*x^3 + 8/3*x^2 - 1/6*x,  45/4*x^6 - 75/2*x^5
+ 45*x^4 -137/6*x^3 + 17/4*x^2 - 1/7*x], (x,0,1))
show(P, frame=True, title=
'Polynomials with p_n(0) = 0 and p_n(1) = B(n) = BernoulliNumber(n)')

Other sequences include: A007840 associated with the ordered Stirling cycle numbers A225479 and generated by "f(n) = (n-1)/n if n>1 else -1" and the Fubini numbers A000670 associated with the ordered Stirling set numbers A131689 (A019538) and generated by "f(n) = 1/n if n > 1 else -1". The reader might also consider the related pair A006153 (sequence) and A199673 (triangle) generated by "f(n) = 1/(n-1) if n > 1 else -1".

The clause "if n > 1 else -1" in the above examples actual means that we are evaluating the polynomials at x = -1 instead of x = 1. This is a simple but noteworthy trick which can be described in full generality as: changing the sign of f(1) turns the row sums into the alternating row sums.

♦ Recurrences for the general case

Since the number of partitions increase rapidly a recurrence for the partition transformation is a necessary condition for an efficient calculation when n gets large. Luckily partitions obey beautiful recurrences and basically we only have to translate these recurrences to our setup.

def PartialPtransRec(n, f, norm = None):
    i = 1; F = [1] 
    while i <= n: F.append(F[i-1]*f(i)); i += 1
    @cached_function
    def PRec(n, k):
        if k==0 and n==0: return 1
        if k==0: return 0
        if k==n: return -PRec(n-1,n-1)*F[1]
        return -sum(F[i]*PRec(n-i,k-1) for i in (1..n-k+1))
    R = [PRec(n, k) for k in (0..n)]
    if norm == None: return R
    return [norm(n,k)*R[k] for k in (0..n)]

For example the Lah numbers of order 2 which we will discuss below can be computed as:

lah2 = lambda n: 1 if n == 1 else ((n-1)^2+1)/(n*(4*n-2)) 
norm = lambda n,k: (-1)^k*factorial(2*n)/factorial(2*k)
for n in (0..6): 
    print PartialPtransRec(n, lah2, norm)

[1]
[0, 1]
[0, 2, 1]
[0, 10, 10, 1]
[0, 100, 140, 28, 1]
[0, 1700, 2900, 840, 60, 1]
[0, 44200, 85800, 31460, 3300, 110, 1]

We can now compute also the inverse P-transformation efficiently by applying inversion to this recursive variant of the P-transformation. But we can do even better: we can also calculate the inverse of the P-transformation directly by a recurrence.

def InversePartialPtransRec(n, f, norm = None):
    i = 1; F = [1] 
    while i <= n: F.append(F[i-1]*f(i)); i += 1
    @cached_function
    def PRec(n, k):
        if k==0 and n==0: return 1
        if k==0: return 0
        if k==n: return -PRec(n-1,n-1)/F[1]
        return -(PRec(n-1,k-1)+sum(F[i]*PRec(n,k+i-1) for i in (2..n-k+1)))/F[1]
    R = [PRec(n, k) for k in (0..n)]
    if norm == None: return R
    return [norm(n,k)*R[k] for k in (0..n)]

As an example we invent an interesting triangle not in the OEIS related to the Euler numbers and which essentially is the inverse of Joffe's central differences of zero:

eul = lambda n: 1/((2*n-1)*(2*n))
nrm = lambda n,k: factorial(2*n)/4^k

for n in (0..6):
    print InversePartialPtransRec(n, eul, nrm)
    
[1]
[0, -1]
[0, -2, 6]
[0, -16, 60, -90]
[0, -288, 1176, -2520, 2520]
[0, -9216, 39360, -98280, 151200, -113400]    

The significance of this triangle is unexplored, though.

♦ Benchmark

Let's check whether our effort has paid off.

import time

def timeit(statement):
    durations = []
    for i in range(3):
        start = time.time()
        eval(statement)
        end = time.time()
        t = end - start
        durations.append(t)
    return min(durations)

def benchmark(stmt):
    t = timeit(stmt)
    r = round(t*1000, 3)
    print "%s best time: %s ms" % (stmt, r)

for n in range(10,50,10):
    benchmark("PartialPtrans(%d, lambda n: n)" % n)

PartialPtrans(10, lambda n: n) best time: 28.879 ms
PartialPtrans(20, lambda n: n) best time: 591.164 ms
PartialPtrans(30, lambda n: n) best time: 7204.341 ms
PartialPtrans(40, lambda n: n) best time: 57720.778 ms

for n in range(10,50,10):
    benchmark("PartialPtransRec(%d, lambda n: n)" % n)

PartialPtransRec(10, lambda n: n) best time: 10.783 ms
PartialPtransRec(20, lambda n: n) best time: 33.621 ms
PartialPtransRec(30, lambda n: n) best time: 80.299 ms
PartialPtransRec(40, lambda n: n) best time: 161.223 ms

Although the algorithm now looks quite efficient the data structure used is not satisfactory. In the typical use case one wants to compute the first n rows of the lower triangular matrix at once, and not row by row. Therefore we give yet another function which does this and additionally integrates the inverse P-transformation and the reduced P-transformation (which we will discuss in the next section) in a single function.

def PtransMatrix(dim, f, norm = None, inverse = False, reduced = False):
    i = 1; F = [1] 
    if reduced: 
        while i <= dim: F.append(f(i)); i += 1
    else: 
        while i <= dim: F.append(F[i-1]*f(i)); i += 1
            
    C = [[0 for k in range(m+1)] for m in range(dim)] 
    C[0][0] = 1
    if inverse:
        for m in (1..dim-1): 
            C[m][m] = -C[m-1][m-1]/F[1]
            for k in range(m-1, 0, -1):     
                C[m][k] = -(C[m-1][k-1]+sum(F[i]*C[m][k+i-1] 
                          for i in (2..m-k+1)))/F[1]  
    else:
        for m in (1..dim-1): 
            C[m][m] = -C[m-1][m-1]*F[1]
            for k in range(m-1, 0, -1):     
                C[m][k] = -sum(F[i]*C[m-i][k-1] for i in (1..m-k+1))
                
    if norm == None: return C 
    for m in (1..dim-1):
        for k in (1..m): C[m][k] *= norm(m,k)
    return C

Even the performance got better:

for n in range(10,50,10):
    benchmark("PtransMatrix(%d, lambda n: n)" % n)
    
PtransMatrix(10, lambda n: n) best time: 6.345 ms
PtransMatrix(20, lambda n: n) best time: 31.637 ms
PtransMatrix(30, lambda n: n) best time: 74.013 ms
PtransMatrix(40, lambda n: n) best time: 140.953 ms    

♦ The reduced P-transformation

The recurrences as implemented not only constitute an efficient mean to compute the P-transformation, they also provide new insight in the structure of the P-transformation. To that end consider the first two lines in the implementation of the recurrences above:

     i = 1; F = [1]; while i <= n: F.append(F[i-1]*f(i)); i += 1

This effectively first transforms the given sequence 'f' into its partial product 'F' before the recursion starts. Thus the P-transformation can be written as

${\mathcal {P}}_{n}(f)={\tilde {\mathcal {P}}}_{n}\left(\prod _{j=1}^{n}f_{j}\right).$

We call ${\tilde {\mathcal {P}}}$ the reduced P-transformation. It starts like:

Obviously this leads to a massive reduction of complexity compared to the original form of the P-transformation.

In this post we will not explore the reduced form although it holds some surprises. For example the identity n → n is mapped by the inverse reduced transform to Catalan's triangle A039598 which counts by a comment of Abdullahi Umar the number of order-preserving transformations of an n-chain with exactly k fixed points.

In any case one should remember that the P-transformations puts four arrows in the quiver to calculate sequences efficiently: the actual P-transformation, its inverse, and the two corresponding reduced transformations.

♦ The coefficients of the multivariate P-polynomials

Let's come back briefly to the multivariate partial P-polynomials. This time we want to extract the coefficients of these polynomials.

def PMultiCoefficients(dim, norm = None, inverse = False):
    def coefficient(p):  # admittedly somewhat obscure 
        c = SR(p).fraction(ZZ).numerator().coefficients() 
        return [0] if not c else c

    f = lambda n: var('x'+str(n))
    P = PtransMatrix(dim, f, norm, inverse)
    return [[coefficient(p) for p in L] for L in P] 

In the ordinary case we get:

PMultiCoefficients(8)
[[1]],
[[0], [-1]],
[[0], [-1], [1]],
[[0], [-1], [2], [-1]],
[[0], [-1], [1, 2], [-3], [1]],
[[0], [-1], [2, 2], [-3, -3], [4], [-1]],
[[0], [-1], [1, 2, 2], [-1, -6, -3], [6, 4], [-5], [1]],
[[0], [-1], [2, 2, 2], [-3, -3, -6, -3], [4, 12, 4], [-10, -5], [6], [-1]]

And in the case of the inverse multivariate partial P-polynomials:

PMultiCoefficients(7, inverse = True)
[[1]],
[[0], [-1]],
[[0], [-1], [1]],
[[0], [-2, 1], [2], [-1]],
[[0], [-5, 5, -1], [5, -2], [-3], [1]],
[[0], [-14, 21, -3, -6, 1], [14, -12, 2], [-9, 3], [4], [-1]]],
[[0], [-42,84,-28,-28,7,7,-1],[42,-56,7,14,-2],[-28,21,-3],[14,-4],[-5],[1]]

Both triangles were not in the OEIS, but column 1 of the inverse P-matrix is Wolfdieter Lang's A111785 named "Array used for power series inversion" (with a different order). The situation is similar to the case of the inverse Bell-transform where previously also only column 1 could be found in the database but not the entire matrix of the inverse coefficients. (The triangles are now A269941 and A269942.)

But what happens if we replace the sublists by their sums? We get the triangle below, and in fact we get it in both cases!

[[1],
 [0, -1],
 [0, -1, 1],
 [0, -1, 2, -1],
 [0, -1, 3, -3, 1],
 [0, -1, 4, -6, 4, -1],
 [0, -1, 5, -10, 10, -5, 1],
 [0, -1, 6, -15, 20, -15, 6, -1],
 [0, -1, 7, -21, 35, -35, 21, -7, 1]]

♦  Prototypic examples

As prototypical examples for the P transformation we choose the unsigned Lah numbers, the Stirling set and the (unsigned) Stirling cycle numbers and Kronecker's delta. We show these examples in the form taken by the partial P-polynomials. This time we use a more mathematical notation. So let's define P^k_n(f) as the P transform of f restricted to the partitions of n with largest part k; this is exactly what the partial P-polynomials generated by f represent. Then:

Signless Lah numbers

$(-1)^{k}{\frac {n!}{k!}}\,{\mathcal {P}}_{n}^{k}\left(1,1,1,\ldots \right)=\left|{n \atop k}\right|$

Stirling cycle numbers

$(-1)^{k}{\frac {n!}{k!}}\,{\mathcal {P}}_{n}^{k}\left(1,{\frac {1}{2}},{\frac {2}{3}},\ldots \right)=\left[{n \atop k}\right]$

Stirling set numbers

$(-1)^{k}{\frac {n!}{k!}}\,{\mathcal {P}}_{n}^{k}\left(1,{\frac {1}{2}},{\frac {1}{3}},\ldots \right)=\left\{{n \atop k}\right\}$

Kronecker delta

$(-1)^{k}{\frac {n!}{k!}}\,{\mathcal {P}}_{n}^{k}(1,0,0,\ldots )=\delta _{n,k}$

There is a natural way to generalize the unsigned Lah numbers, the Stirling cycle numbers and the Stirling set numbers to a full hierarchy of combinatorial numbers in a uniform way (I described this elsewhere). The case m = 2 in this hierarchy gives for the Stirling set numbers the central factorial numbers T(2n, 2k) in Riordan's notation (A036969) and for the Stirling cycle numbers the central factorial numbers t(2n, 2k) (A204579). The case of the Lah numbers has not been considered before as far as I know (the numbers are now registered as A268434). All these cases can be computed with the partial P-polynomials as easily as the base cases.

The generators and the normalization:

stirset2   = lambda n: 1 if n == 1 else           1/(n*(4*n-2)) 
stircycle2 = lambda n: 1 if n == 1 else     (n-1)^2/(n*(4*n-2)) 
lah2       = lambda n: 1 if n == 1 else ((n-1)^2+1)/(n*(4*n-2)) 
norm       = lambda n,k: (-1)^k*factorial(2*n)/factorial(2*k)

Stirling set order 2

M = PtransMatrix(7, stirset2, norm) 
for m in M: print m    
[1]
[0, 1]
[0, 1, 1]
[0, 1, 5, 1]
[0, 1, 21, 14, 1]
[0, 1, 85, 147, 30, 1]
[0, 1, 341, 1408, 627, 55, 1]

M = PtransMatrix(7, stirset2, norm, inverse = True) 
for m in M: print m     
[1]
[0, 1]
[0, 1, 1]
[0, 4, 5, 1]
[0, 36, 49, 14, 1]
[0, 576, 820, 273, 30, 1]
[0, 14400, 21076, 7645, 1023, 55, 1]

Stirling cycle order 2

M = PtransMatrix(7, stircycle2, norm)   # A204579
for m in M: print m     
[1],
[0, 1],
[0, 1, 1],
[0, 4, 5, 1],
[0, 36, 49, 14, 1],
[0, 576, 820, 273, 30, 1],
[0, 14400, 21076, 7645, 1023, 55, 1]  
 
M = PtransMatrix(7, stircycle2, norm, inverse = True) 
for m in M: print m 
[1]
[0, 1]
[0, 1, 1]
[0, 1, 5, 1]
[0, 1, 21, 14, 1]
[0, 1, 85, 147, 30, 1]
[0, 1, 341, 1408, 627, 55, 1] 

Lah order 2

M = PtransMatrix(7, lah2, norm) 
for m in M: print m 
[1]
[0,     1]
[0,     2,     1]
[0,    10,    10,     1]
[0,   100,   140,    28,    1]
[0,  1700,  2900,   840,   60,   1]
[0, 44200, 85800, 31460, 3300, 110, 1]

The reader might note that the triangles stirset2 and stircycle2 are inverse to each other and the lah2 triangle is self-inverse (as it should!).

Signless Lah numbers of order 2

$(-1)^{k}{\frac {(2n)!}{(2k)!}}\,{\mathcal {P}}_{n}^{k}\left(1,{\frac {1}{6}},{\frac {1}{6}},{\frac {5}{28}},\ldots \right)=\left|{n \atop k}\right|^{[2]}$

Stirling cycle numbers of order 2 (a.k.a. central factorial numbers t(2n, 2k))

$(-1)^{k}{\frac {(2n)!}{(2k)!}}\,{\mathcal {P}}_{n}^{k}\left(1,{\frac {1}{12}},{\frac {2}{15}},{\frac {9}{56}},\ldots \right)=\left[{n \atop k}\right]^{[2]}$

Stirling set numbers of order 2 (a.k.a. central factorial numbers T(2n, 2k))

$(-1)^{k}{\frac {(2n)!}{(2k)!}}\,{\mathcal {P}}_{n}^{k}\left(1,{\frac {1}{12}},{\frac {1}{30}},{\frac {1}{56}},\ldots \right)=\left\{{n \atop k}\right\}^{[2]}$

But why call $\textstyle \left|{n \atop k}\right|^{[2]}$ the Lah numbers of order 2? Well, there are several convincing arguments that this is the 'right' generalization of the Lah numbers. We indicate just the simplest one: Lah numbers of order 1 tie in with the Stirling numbers of order 1 by

$\left|{n \atop k}\right|=\sum \limits _{j\,=\,k}^{n}\left[{n \atop j}\right]\left\{{j \atop k}\right\}\,,$

and Lah numbers of order 2 have the representation

$\left|{n \atop k}\right|^{[2]}=\sum \limits _{j\,=\,k}^{n}\left[{n \atop j}\right]^{[2]}\left\{{j \atop k}\right\}^{[2]}\,.$

♦  Stirling and Lah numbers of higher order

We are inclined to say that a generalization of the Stirling numbers is only satisfying if it leads simultaneously to a satisfactory generalization of the Lah numbers. To see the full picture, let us consider briefly the third step in the hierarchy, but this time focusing on the structure of the recurrence relation behind these numbers.

def T(n, k, w):
    if n==k: return 1 
    if k<0 or k>n: return 0
    return T(n-1, k-1, w) + w(n,k)*T(n-1, k, w)

Choosing for w the weight functions stirset3, stircycle3 and lah3 gives us the desired generalizations.

stirset3   = lambda n,k: k^3
stircycle3 = lambda n,k: (n-1)^3
lah3       = lambda n,k: stircycle3(n,k)+stirset3(n,k)

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

[1]
[0, 1]
[0, 2, 1]
[0, 18, 18, 1]
[0, 504, 648, 72, 1]
[0, 32760, 47160, 7200, 200, 1]
[0, 4127760, 6305040, 1141560, 45000, 450, 1]
[0, 895723920, 1416456720, 283704120, 13741560, 198450, 882, 1]

Replacing the '3's in the weight functions by a nonnegative m now gives the full generalization. In particular we see that the case m = 0 reduces both kinds of Stirling numbers to Pascal's triangle and that the Lah numbers of order 0 are given by A038207, which is the square of Pascal's triangle.

The hierarchy of the Stirling/Lah numbers.

Order	Stirling set	Stirling cycle	Lah numbers
0	A007318	A007318	A038207
1	A048993	A132393	A271703
2	A269945	A269944	A268434
3	A269948	A269947	A269946

♦  Some additional Stirling number identities

We now use the P-transformation to introduce two number triangles which give raise to six apparently new identities related to the Stirling numbers.

Stirling set^* numbers (fictitious name)

$\left\{{n \atop k}\right\}^{\star }=(-1)^{k}(2n)!\,{\mathcal {P}}_{n}^{k}\left(1,{\frac {1}{2}},{\frac {1}{3}},{\frac {1}{4}},\ldots \right)$

Stirling cycle^* numbers (fictitious name)

$\left[{n \atop k}\right]^{\star }=(-1)^{k}(2n)!\,{\mathcal {P}}_{n}^{k}\left(1,{\frac {1}{2}},{\frac {2}{3}},{\frac {3}{4}},\ldots \right)$

The numerical values of these numbers can be found in A268437 and A268438. The identities display a nice formal symmetry:

The values are in A268435 and A268436. In order to obtain another representation of these numbers we look at two more number triangles:

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

$\operatorname {V} _{n,k}\,=\,\sum _{m=0}^{k}(-1)^{m+k}{\binom {n+k}{n+m}}\,\left[{n+m \atop m}\right]$

With this we get

$\left\{{n \atop k}\right\}^{\star }={\frac {(2n)!}{(n+k)^{\underline {n}}}}\,\operatorname {W} _{n,k},\qquad \left[{n \atop k}\right]^{\star }={\frac {(2n)!}{(n+k)^{\underline {n}}}}\,\operatorname {V} _{n,k}.$

The numbers W_n,k are known as Ward numbers. Charles Jordan gives in his 'Calculus of Finite Differences', Chelsea 1950, short tables of these numbers on page 152 and page 172 respectively. Compare A269939 and A269940.

♦  Maple code

First we want to check the formulas in the last section for some small values.

C := (n,k) -> binomial(n,k):
R := (n,k) -> pochhammer(n,k):
F := (n,k) -> (-1)^k*R(-n,k):
N := (n,k) -> ((2*n)!/F(n+k,n)):
S := (n,k) -> Stirling2(n,k):
s := (n,k) -> abs(Stirling1(n,k)):

B := (n,k) -> N(n,k)*add((-1)^(m+k)*C(n+k,n+m)*S(n+m,m), m=0..k):
b := (n,k) -> N(n,k)*add((-1)^(m+k)*C(n+k,n+m)*s(n+m,m), m=0..k):

A1 := (n,k) -> R(n-k+1,n-k)*S(n,k): 
A2 := (n,k) -> C(n,k)*add(C(k,i)*B(n-k,i),i=0..k):
A3 := (n,k) -> C(-k,-n)*add(C(-n,i)*b(n-k,i),i=0..n-k): 
  
B1 := (n,k) -> R(n-k+1,n-k)*s(n,k): 
B2 := (n,k) -> C(n,k)*add(C(k,i)*b(n-k,i), i=0..k):
B3 := (n,k) -> C(-k,-n)*add(C(-n,i)*B(n-k,i), i=0..n-k):

Here the numbers 'B' and 'b' are the stared Stirling numbers using the given representations W_n,k and V_n,k and A1, A2, A3 the expressions in the first line of identities and B1, B2, B3 the expressions in the second line of identities.

test := f -> seq(print(seq(f(n,k),k=0..n)),n=0..6):
test(A1);test(A2);test(A3);
test(B1);test(B2);test(B3);

Our test shows that A1, A2 and A3 are equal to A268435 and B1, B2 and B3 equal to A268436 in the tested range as expected.

It is a disturbing fact that SageMath has only implemented a crippled form of the binomial coefficients which does not perform this test correctly. So we have to say that SageMath is not a suitable tool for the calculation of many sequences in OEIS -- a sad finding.

Now let's come back to the partial P-transformation and the P-transformation:

partial_ptrans := proc(n, f:=NULL) local q,p,r,g,R; 
if n = 0 then return [1] fi;
if f = NULL then g := n -> x[n] else g := f fi;
R := [seq(0,j=0..n)]; 
for q in combinat:-partition(n) do
   p := [op(ListTools:-Reverse(q)),0]; r := p[1]+1;
   mul(binomial(p[j], p[j+1])*g(j)^p[j], j=1..nops(q));
   R[r] := R[r] - (-1)^r*%; 
od; R end:

p_trans := proc(n,f:=NULL) add(k,k=partial_ptrans(n,f)) end:

Example calls:

partial_ptrans(5); 
partial_ptrans(6, m->m);
seq(p_trans(n, m->1/(m+1))*n!, n=0..10);

[0, -x[1] x[2] x[3] x[4] x[5], 
 2x[1]^2 x[2] x[3] x[4] + 2x[1]^2 x[2]^2 x[3], 
-3x[1]^3 x[2] x[3] -3x[1]^3 x[2]^2, 4x[1]^4 x[2], -x[1]^5 ]

[0, -720, 372, -152, 48, -10, 1]

1, -1/2, 1/6, 0, -1/30, 0, 1/42, 0, -1/30, 0, 5/66

The P-matrix with parameter dimension dim and the optional parameters generator and normalisator.

p_matrix := proc(dim, generator:=NULL, normal:=NULL) 
local n, k, w, M, T; 
M := Matrix(dim, shape=triangular[lower]); M[1,1] := 1;
if normal=NULL then w := (n,k)->1 else w:=normal fi;
for n from 1 to dim do 
    T := partial_ptrans(n-1, generator);
    for k from 1 to n do M[n,k] := w(n-1,k-1)*T[k] od; 
od; M end:

Example calls:

p_matrix(6);
p_matrix(6, m->m);

f := n -> `if`(n=1,1,1/(n*(4*n-2))):
w := (n,k) -> (-1)^k*(2*n)!/(2*k)!:
p_matrix(6, f, w);

The sum of coefficients of the generated polynomials:

p_coeffsum := proc(len, f) local R,C,n,k;
R := Array(0..len-1); C := Array(0..len-1); 
R[0] := 1; C[0] := 1; 
for n from 1 to len-1 do
    for k from n by -1 to 1 do C[k] := C[k-1]*f(k) od;
    C[0] := -add(C[k], k=1..n); R[n] := %; od; 
convert(R,list) end:

p_coeffsum(9, m->m);
[1, -1, -1, -3, -13, -71, -461, -3447, -29093]

♦  Summary

We have introduced the P-transformation, a sequence-to-triangle transformation, based on the partitions of n with largest part k. This transformation is — next to the closely related Bell-transformation — one of the core transformations defined on sequences. That this transformation has not received the attention it deserves in the transformation-libraries which are distributed with the OEIS is in our opinion due to the fact that important sequences are generated through rational sequences. The Euler and the Bernoulli numbers are outstanding examples for this point of view and prove the usefulness of this transformation.

♦ Some P-polynomials evaluated at x=1 or x=-1

A000007	A000670	A001896	A003319	A006153	A006232
A006252	A006568	A006569	A007840	A027641	A027642
A036280	A036281	A036968	A051296	A075178	A077607
A089148	A101686	A113871	A118196	A135920	A154288
A154289	A167894	A198631	A226158	A248964	A249024

---

Accompanying, there is a SageMath Jupyter Notebook which can be downloaded from GitHub.

Thanks to my editor Hans Havermann!