This site is supported by donations to The OEIS Foundation.
User:Peter Luschny/PTransform
The PTransform
Contents
 1 ♦ Ppolynomials
 2 ♦ Representations of the Euler and Bernoulli numbers
 3 ♦ The inverse of the Ptransform
 4 ♦ Recurrence for the case x = 1
 5 ♦ Recurrences for the general case
 6 ♦ Benchmark
 7 ♦ The reduced Ptransformation
 8 ♦ The coefficients of the multivariate Ppolynomials
 9 ♦ Prototypic examples
 10 ♦ Stirling and Lah numbers of higher order
 11 ♦ Some additional Stirling number identities
 12 ♦ Maple code
 13 ♦ Summary
 14 ♦ Some Ppolynomials evaluated at x=1 or x=1
♦ Ppolynomials
We introduce a sequencetotriangle transformation which we will call the partitiontransformation (or short Ptransformation) 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 Ptransformation 'PartialPtrans' and the Ptransformation '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_{1}, x_{2}, x_{3}, ....
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 Ppolynomials' to distinguish them from the partial Ppolynomials, 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 Ppolynomial is just the sum of the corresponding partial Ppolynomials 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_{1}, x_{2}, x_{3}, ... by the values f(1), f(2), f(3), ... (in other words, it evaluates the (multivariate) Ppolynomials 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(Len2): print AsPolynomials(n, f, norm) print "Coefficients:" for n in range(Len2): 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 Ptransformation. 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 Ptransformations 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 Ptransformation is interesting we show how the Euler and the Bernoulli numbers are represented as Ptransforms. These examples also show that when we are speaking of a sequencetotriangle 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/((2n1)(2n)) and normalize the result by multiplying with (2n)!. Then we evaluate.
f = lambda n: 1/((2*n1)*(2*n)) norm = lambda n,k: factorial(2*n) Len = 8 print "As polynomials:" for n in range(Len2): print AsPolynomials(n, f, norm) print "Coefficients:" for n in range(Len2): 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/(2cosh(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 Ptransformation 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:

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 Ppolynomials 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)!/(22^(2n)). Then we evaluate:
f = lambda n: 1/((2*n)*(2*n+1)) norm = lambda n,k: factorial(2*n)/(22^(2*n))
print "\n # As polynomials:" for n in range(Len2): print AsPolynomials(n, f, norm) print "\n # Coefficients" for n in range(Len2): 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 nonvanishing constant terms for even n.
The finestructure of the Bernoulli numbers which we used in this calculation written in mathematical notation:

♦ The inverse of the Ptransform
The inverse of the Ptransform is defined as the inverse of the lower triangular matrix of the partial Ppolynomials.
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..dim1): M[n][n] = 1/A[n][n] for k in range(n1, 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..dim1): for k in (1..n): M[n][k] *= norm(n,k) return M
In its full generality the matrix of the inverse partial Ppolynomials 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^4x2*x3*x4/x1^4, 5*x2^2/x1^42*x2*x3/x1^4, 3*x2/x1^4, 1/x1^4]
We see that x_{1} != 0 is a necessary condition for the existence of the inverse Ptransform. Column 1 of this matrix will be called the Pinverse of the sequence x_{1}, x_{2}, x_{3},... . The inverse Pmatrix can be cleaned up a bit if we assume x_{1} = 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 Pinverse 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..n1)] print InversePsequence(8, lambda n: n) [1, 2, 2, 4, 4, 48, 336]
♦ Recurrence for the case x = 1
Evaluating the Ppolynomials 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]*(len1) for n in (1..len1): for k in range(n, 0, 1): C[k] = C[k1] * 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]*(len1) for n in (1..len1): for k in range(n, 0, 1): C[k] = C[k1] / (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]*(len1) for n in (1..len1): f *= n for k in range(n, 0, 1): C[k] = C[k1] / (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) = (n1)/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/(n1) 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[i1]*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(n1,n1)*F[1] return sum(F[i]*PRec(ni,k1) for i in (1..nk+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 ((n1)^2+1)/(n*(4*n2)) 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 Ptransformation efficiently by applying inversion to this recursive variant of the Ptransformation. But we can do even better: we can also calculate the inverse of the Ptransformation directly by a recurrence.
def InversePartialPtransRec(n, f, norm = None): i = 1; F = [1] while i <= n: F.append(F[i1]*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(n1,n1)/F[1] return (PRec(n1,k1)+sum(F[i]*PRec(n,k+i1) for i in (2..nk+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*n1)*(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 Ptransformation and the reduced Ptransformation (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[i1]*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..dim1): C[m][m] = C[m1][m1]/F[1] for k in range(m1, 0, 1): C[m][k] = (C[m1][k1]+sum(F[i]*C[m][k+i1] for i in (2..mk+1)))/F[1] else: for m in (1..dim1): C[m][m] = C[m1][m1]*F[1] for k in range(m1, 0, 1): C[m][k] = sum(F[i]*C[mi][k1] for i in (1..mk+1)) if norm == None: return C for m in (1..dim1): 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 Ptransformation
The recurrences as implemented not only constitute an efficient mean to compute the Ptransformation, they also provide new insight in the structure of the Ptransformation. 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[i1]*f(i)); i += 1
This effectively first transforms the given sequence 'f' into its partial product 'F' before the recursion starts. Thus the Ptransformation can be written as
We call the reduced Ptransformation. It starts like:
Obviously this leads to a massive reduction of complexity compared to the original form of the Ptransformation.
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 orderpreserving transformations of an nchain with exactly k fixed points.
In any case one should remember that the Ptransformations puts four arrows in the quiver to calculate sequences efficiently: the actual Ptransformation, its inverse, and the two corresponding reduced transformations.
♦ The coefficients of the multivariate Ppolynomials
Let's come back briefly to the multivariate partial Ppolynomials. 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 Ppolynomials:
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 Pmatrix 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 Belltransform 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 Ppolynomials. 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 Ppolynomials generated by f represent. Then:
Signless Lah numbers
Stirling cycle numbers
Stirling set numbers
Kronecker delta
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 Ppolynomials as easily as the base cases.
The generators and the normalization:
stirset2 = lambda n: 1 if n == 1 else 1/(n*(4*n2)) stircycle2 = lambda n: 1 if n == 1 else (n1)^2/(n*(4*n2)) lah2 = lambda n: 1 if n == 1 else ((n1)^2+1)/(n*(4*n2)) 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 selfinverse (as it should!).
Signless Lah numbers of order 2
Stirling cycle numbers of order 2 (a.k.a. central factorial numbers t(2n, 2k))
Stirling set numbers of order 2 (a.k.a. central factorial numbers T(2n, 2k))
But why call 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
and Lah numbers of order 2 have the representation
♦ 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(n1, k1, w) + w(n,k)*T(n1, 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: (n1)^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.
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 Ptransformation to introduce two number triangles which give raise to six apparently new identities related to the Stirling numbers.
Stirling set^{*} numbers (fictitious name)
Stirling cycle^{*} numbers (fictitious name)
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:
With this we get
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(nk+1,nk)*S(n,k): A2 := (n,k) > C(n,k)*add(C(k,i)*B(nk,i),i=0..k): A3 := (n,k) > C(k,n)*add(C(n,i)*b(nk,i),i=0..nk): B1 := (n,k) > R(nk+1,nk)*s(n,k): B2 := (n,k) > C(n,k)*add(C(k,i)*b(nk,i), i=0..k): B3 := (n,k) > C(k,n)*add(C(n,i)*B(nk,i), i=0..nk):
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 Ptransformation and the Ptransformation:
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 Pmatrix 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(n1, generator); for k from 1 to n do M[n,k] := w(n1,k1)*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*n2))): 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..len1); C := Array(0..len1); R[0] := 1; C[0] := 1; for n from 1 to len1 do for k from n by 1 to 1 do C[k] := C[k1]*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 Ptransformation, a sequencetotriangle transformation, based on the partitions of n with largest part k. This transformation is — next to the closely related Belltransformation — one of the core transformations defined on sequences. That this transformation has not received the attention it deserves in the transformationlibraries 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 Ppolynomials evaluated at x=1 or x=1
Accompanying, there is a Jupyter Notebook which was developed on SageMathCloud. SageMath is a free open source alternative to Magma, Maple, Mathematica and Matlab. Get a free account on SageMathCloud (SMC), download this notebook, upload it to SMC and start exploring!
Thanks to my editor Hans Havermann!