This site is supported by donations to The OEIS Foundation.
User:Peter Luschny/BellTransform
The Bell Transform
Summary: We introduce the partial Bell polynomials in two flavors: in the traditional way based on the variables x1, x2, ... and based on the variables x0, x1, x2, .... The coefficients of the x0-based univariate polynomials give the new and unexplored A257563.
Next we base the computation of the Bell transform on partitions which gives a surprising fast algorithm for small n which is similar to the algorithm used in SageMath for Bell polynomials.
Iterating the Bell transform using the row sums of the triangles as input for the next level gives us the higher order Bell numbers. In particular we recover Joerg Arndt's award winning A187761 as the second-order Bell numbers.
We then apply the Bell transform and its inverse to the power functions, the multifactorials, the rising factorials and the falling factorials showing more than 50 examples of sequences in the database. We found that almost none of them was described as the result of a Bell transformation.
Contents
- 1 ♦ Bell polynomials
- 2 ♦ The Bell transform based on partitions
- 3 ♦ The Bell matrix
- 4 ♦ The inverse Bell transform
- 5 ♦ The Bell inverse of a sequence
- 6 ♦ The coefficients of the Bell polynomials
- 7 ♦ Bell numbers of higher order
- 8 ♦ Associated Stirling subset numbers
- 9 ♦ The Bell transform of multifactorials
- 10 ♦ The Bell transform of powers
- 11 ♦ The Bell transform of rising factorials
- 12 ♦ The Bell transform of falling factorials
- 13 ♦ The Bell transform of monotonic factorials
- 14 ♦ An implementation for Maple
- 15 ♦ Index of triangles generated by the Bell transform
- 16 ♦ References
♦ Bell polynomials
The Bell transform is a transformation which maps a sequence
a = a_0, a_1, a_2,...
to a triangle
T(0,0)
T(1,0) T(1,1)
T(2,0) T(2,1) T(2,2)
T(3,0) T(3,1) T(3,2) T(3,3)
defined by
T(n,k) = Sum_{j=0..n-k+1} binomial(n-1,j-1)*a_j*T(n-j,k-1)
if k != 0 else T(n, k) = (a_0)^n.
With Sage this can be implemented as
def bell_trans(n,k,a):
@cached_function
def T(n, k):
if k==0: return a[0]^n
return sum(binomial(n-1,j-1)*a[j]*T(n-j,k-1) for j in (0..n-k+1))
return T(n,k)
Let us look at the first few rows of the generated triangle.
(*)
a_0^0
a_0^1, a_1
a_0^2, a_1*a_0+a_2, a_1^2
a_0^3, a_1*a_0^2+2*a_2*a_0+a_3, a_1^2*a_0+3*a_2*a_1, a_1^3
a_0^4, a_1*a_0^3+3*a_2*a_0^2+3*a_3*a_0+a_4, a_1^2*a_0^2+5*a_1*a_2*a_0
+4*a_3*a_1+3*a_2^2, a_1^3*a_0+6*a_2*a_1^2, a_1^4
For example to compute the unsigned Stirling numbers of the first kind A132393 one can write
s = [0]+[factorial(i) for i in (0..5)]
for n in (0..5): [bell_trans(n,k,s) for k in (0..n)]
[1]
[0, 1]
[0, 1, 1]
[0, 2, 3, 1]
[0, 6, 11, 6, 1]
[0, 24, 50, 35, 10, 1]
A second example shows that there are other important numbers which can be generated this way which do not always reduce to integers.
N = 6
t = [2^(1-j) if is_odd(j) else 0 for j in (0..N)]
for n in (0..N): [bell_trans(n,k,t) for k in (0..n)]
[1]
[0, 1]
[0, 0, 1]
[0, 1/4, 0, 1]
[0, 0, 1, 0, 1]
[0, 1/16, 0, 5/2, 0, 1]
[0, 0, 1, 0, 5, 0, 1]
These are the coefficients of the central factorials [1]
Another way to look at the symbolic sums in the triangle (*) above is to interpret the an as variables xn of a multivariate polynomial. Then one arrives at the partial Bell polynomials:
def partial_bell_polynomial(n,k):
v = var(['x_'+str(i) for i in (0..n+1)])
return bell_trans(n,k,v).expand()
for n in (0..4):
[partial_bell_polynomial(n,k) for k in (0..n)]
[1]
[x_0, x_1]
[x_0^2, x_0*x_1 + x_2, x_1^2]
[x_0^3, x_0^2*x_1 + 2*x_0*x_2 + x_3, x_0*x_1^2 + 3*x_1*x_2, x_1^3]
[x_0^4, x_0^3*x_1 + 3*x_0^2*x_2 + 3*x_0*x_3 + x_4, x_0^2*x_1^2 +
5*x_0*x_1*x_2 + 3*x_2^2 + 4*x_1*x_3, x_0*x_1^3 + 6*x_1^2*x_2, x_1^4]
Note that these polynomials depend on the variables x0, x1, x2, ... For some but probably not for a good reason the partial Bell polynomials are often defined to depend only on the variables x1, x2, ... In this case the first column of the above triangle is just ignored (except for the case (0, 0) which is added somewhat unmotivated and unsystematic). To get this in the older literature common form (see for example L. Comtet, Advanced Combinatorics, page 135) from our form is easy: we just have to set x0 = 0.
With Sage we can write:
for n in (0..4):
[partial_bell_polynomial(n,k).subs(x_0=0) for k in (0..n)]
[1]
[0, x_1]
[0, x_2, x_1^2]
[0, x_3, 3*x_1*x_2, x_1^3]
[0, x_4, 3*x_2^2 + 4*x_1*x_3, 6*x_1^2*x_2, x_1^4]
So far we have only looked at the partial Bell polynomials. The Bell polynomials are defined as
B_n = Sum_{k=0..n} B_{n,k}.
def bell_polynomial(n):
return sum(partial_bell_polynomial(n,k) for k in (0..n))
for n in (0..3): bell_polynomial(n)
1
x_0 + x_1
x_0^2 + x_0*x_1 + x_1^2 + x_2
x_0^3 + x_0^2*x_1 + x_0*x_1^2 + x_1^3 + 2*x_0*x_2 + 3*x_1*x_2 + x_3
Two things remain to be mentioned:
First, it is sometimes convenient to reduce the infinite number of variables x_0, x_1, x_2, ... to one variable x only.
Second, in the considerations above we worked with 'symbolic sums' (or in the Sage parlance in the Symbolic Ring SR). But of course we can also work with 'real' polynomials.
Putting these two requirements together we arrive at polynomials in the univariate polynomial ring in x over the rational field.
def univariate_bell_polynomial(n):
p = bell_polynomial(n).subs(x_0=0)
q = p({p.variables()[i]:x for i in range(len(p.variables()))})
R = PolynomialRing(QQ, 'x')
return R(q)
for n in (0..6): univariate_bell_polynomial(n)
1
x
x^2 + x
x^3 + 3*x^2 + x
x^4 + 6*x^3 + 7*x^2 + x
x^5 + 10*x^4 + 25*x^3 + 15*x^2 + x
x^6 + 15*x^5 + 65*x^4 + 90*x^3 + 31*x^2 + x
If we extract the coefficients of these polynomials we find:
for n in (0..6): univariate_bell_polynomial(n).list()
[1]
[0, 1]
[0, 1, 1]
[0, 1, 3, 1]
[0, 1, 7, 6, 1]
[0, 1, 15, 25, 10, 1]
[0, 1, 31, 90, 65, 15, 1]
Lo and behold, these are the Stirling subset numbers A048993! This must have been observed before ;-)
Note the first substitution subs(x_0=0) which indicates the 'classical case' in the sense discussed above. If instead we use the uniform substitution subs(x_n=x) for all n we get:
def x0_based_univariate_bell_polynomial(n):
p = bell_polynomial(n)
q = p({p.variables()[i]:x for i in range(len(p.variables()))})
R = PolynomialRing(QQ,'x')
return R(q)
for n in (0..6): x0_based_univariate_bell_polynomial(n)
1
2*x
3*x^2 + x
4*x^3 + 5*x^2 + x
5*x^4 + 14*x^3 + 10*x^2 + x
6*x^5 + 30*x^4 + 48*x^3 + 19*x^2 + x
7*x^6 + 55*x^5 + 158*x^4 + 149*x^3 + 36*x^2 + x
for n in (0..6): x0_based_univariate_bell_polynomial(n).list()
[1]
[0, 2]
[0, 1, 3]
[0, 1, 5, 4]
[0, 1, 10, 14, 5]
[0, 1, 19, 48, 30, 6]
[0, 1, 36, 149, 158, 55, 7]
The row sums
1, 2, 4, 10, 30, 104, 406, 1754, 8280, 42294, 231950, ...
are Paul Barry's A186021, Bell(n)*(2-0^n). Geoffrey Critzer commented: These are "the number of collections of subsets of {1,2,...,n-1} that are pairwise disjoint."
The coefficients of the x0-based univariate polynomials are now in A257563. What is their combinatorial meaning?
♦ The Bell transform based on partitions
We did not yet consider the question of efficiency of our implementation, to which we turn now.
If we delete the intermediate steps which we introduced for the sake of exposition above we can sum up our implementation as:
def bell_polynomial(n):
X = var(['x_'+str(i) for i in (0..n+1)])
@cached_function
def T(n, k):
if k==0: return k^n
return sum(binomial(n-1,j-1)*T(n-j,k-1)*X[j]
for j in (0..n-k+1)).expand()
return sum(T(n,k) for k in (0..n))
The nice thing is that we can speed up things by explicitly summing over partitions of n instead of using the recursive form. This is somewhat surprising since the relevant formula looks quite formidable; the reader will find it as the first formula on Wikipedia's entry on Bell polynomials.
(Note that Wikipedia calls 'complete Bell polynomials' what we call 'Bell polynomials' and Mathworld calls 'Bell polynomials' what we call 'univariate Bell polynomials'.)
We start with the simpler case of the univariate Bell polynomials and look at an implementation which is essentially the implementation of Sage (GPL licensed), written by Blair Sutton.
def partition_based_univariate_bell_polynomial(n):
polynom = x^n
fn = factorial(n)
for k in (0..n-1):
result = 0
for p in Partitions(n, length=k):
factorial_product = 1
power_factorial_product = 1
for part, count in p.to_exp_dict().iteritems():
factorial_product *= factorial(count)
power_factorial_product *= factorial(part)**count
coefficient = fn//(factorial_product*power_factorial_product)
result += coefficient
polynom += result*x^k
return polynom
So let's benchmark!
%timeit univariate_bell_polynomial(10) 5 loops, best of 3: 105 ms per loop %timeit partition_based_univariate_bell_polynomial(10) 5 loops, best of 3: 14.1 ms per loop %timeit univariate_bell_polynomial(20) 5 loops, best of 3: 632 ms per loop %timeit partition_based_univariate_bell_polynomial(20) 5 loops, best of 3: 79.4 ms per loop %timeit univariate_bell_polynomial(30) 5 loops, best of 3: 2.15 s per loop %timeit partition_based_univariate_bell_polynomial(30) 5 loops, best of 3: 1.16 s per loop
Thus for small n we observe a speedup between 2 and 8.
Back to math. Our main goal is to explore the Bell transform. Therefore we skip the details concerning the implementation of the multivariate Bell polynomials (as polynomials) and give directly the formalism in its most useful form: as a sequence to triangle transformation. To make our setup more generic we now change from lists as input to functions.
def bell_transform(f, n): # partition_based
row = []
fn = factorial(n)
for k in (0..n):
result = 0
for p in Partitions(n, length=k):
factorial_product = 1
power_factorial_product = 1
for part, count in p.to_exp_dict().iteritems():
factorial_product *= factorial(count)
power_factorial_product *= factorial(part)**count
coefficient = fn//(factorial_product*power_factorial_product)
result += coefficient*prod([f(i-1) for i in p])
row.append(result)
return row
♦ The Bell matrix
For example let's look at the family of Stirling-cycle/Lah numbers, which are the Bell transforms of the functions "n → factorial(n+j)" (for some fixed j).
for n in (0..7): bell_transform(factorial, n)
[1]
[0, 1]
[0, 1, 1]
[0, 2, 3, 1]
[0, 6, 11, 6, 1]
[0, 24, 50, 35, 10, 1]
[0, 120, 274, 225, 85, 15, 1]
[0, 720, 1764, 1624, 735, 175, 21, 1]
Unsigned Stirling numbers of the first kind, A132393.
for n in (0..7):
bell_transform(lambda x: factorial(x+1), n)
[1]
[0, 1]
[0, 2, 1]
[0, 6, 6, 1]
[0, 24, 36, 12, 1]
[0, 120, 240, 120, 20, 1]
[0, 720, 1800, 1200, 300, 30, 1]
[0, 5040, 15120, 12600, 4200, 630, 42, 1]
Unsigned Lah numbers, cf. A111596.
for n in (0..7):
bell_transform(lambda x: factorial(x+2), n)
[1]
[0, 2]
[0, 6, 4]
[0, 24, 36, 8]
[0, 120, 300, 144, 16]
[0, 720, 2640, 2040, 480, 32]
[0, 5040, 25200, 27720, 10320, 1440, 64]
[0, 40320, 262080, 383040, 199920, 43680, 4032, 128]
Unsigned coefficients for rewriting generalized falling
factorials into ordinary falling factorials, A136656.
The next step is a trival one: We collect these lists in a matrix, the Bell matrix, and call the function given as input the generator of the Bell matrix. With SageMath we can write this:
def bell_matrix(generator, dim):
row = lambda n: bell_transform(generator, n)
return matrix(ZZ, [row(n)+[0]*(dim-n-1) for n in srange(dim)])
Alternatively, if we want to avoid the matrix constructor and base the Bell transform on its recursive form we can compute the Bell matrix with this function:
def Bell_Matrix(generator, dim):
A = [[0]*(n+1) for n in range(dim)]
for n in range(dim):
A[n][0] = 1 if n==0 else 0
if n>0: A[n][1] = generator(n-1)
for k in (2..n):
A[n][k] = sum(binomial(n-1,j-1)*A[n-j][k-1]*A[j][1] for j in (1..n-k+1))
return A
♦ The inverse Bell transform
The inverse Bell transform is the map generator → inverse(bell_transform(generator)). Here we consider 'bell_transform(generator)' as an infinite lower matrix of which we calculate the matrix inverse.
For example using the factorial function as the generator we can write:
bell_matrix(factorial, 8)
[ 1 0 0 0 0 0 0 0]
[ 0 1 0 0 0 0 0 0]
[ 0 1 1 0 0 0 0 0]
[ 0 2 3 1 0 0 0 0]
[ 0 6 11 6 1 0 0 0]
[ 0 24 50 35 10 1 0 0]
[ 0 120 274 225 85 15 1 0]
[ 0 720 1764 1624 735 175 21 1]
Then the inverse Bell transform can be computed as:
bell_matrix(factorial, 8).inverse()
[ 1 0 0 0 0 0 0 0]
[ 0 1 0 0 0 0 0 0]
[ 0 -1 1 0 0 0 0 0]
[ 0 1 -3 1 0 0 0 0]
[ 0 -1 7 -6 1 0 0 0]
[ 0 1 -15 25 -10 1 0 0]
[ 0 -1 31 -90 65 -15 1 0]
[ 0 1 -63 301 -350 140 -21 1]
Another example are the coefficients of the Abel polynomials A137452 which have the enumeration of the natural numbers as their generator:
bell_matrix(lambda n: n+1, 8).inverse() [ 1 0 0 0 0 0 0 0] [ 0 1 0 0 0 0 0 0] [ 0 -2 1 0 0 0 0 0] [ 0 9 -6 1 0 0 0 0] [ 0 -64 48 -12 1 0 0 0] [ 0 625 -500 150 -20 1 0 0] [ 0 -7776 6480 -2160 360 -30 1 0] [ 0 117649 -100842 36015 -6860 735 -42 1]
This is the inverse of the triangle of idempotent numbers A059297 :
bell_matrix(lambda n: n+1, 8) [ 1 0 0 0 0 0 0 0] [ 0 1 0 0 0 0 0 0] [ 0 2 1 0 0 0 0 0] [ 0 3 6 1 0 0 0 0] [ 0 4 24 12 1 0 0 0] [ 0 5 80 90 20 1 0 0] [ 0 6 240 540 240 30 1 0] [ 0 7 672 2835 2240 525 42 1]
Alternatively, if we want to avoid the matrix constructor and base the inverse Bell transform on its recursive form we can use the following function:
def Inverse_Bell_Matrix(f, dim):
A = Bell_Matrix(f, dim)
M = [[0]*(n+1) for n in range(dim)]
for n in range(dim):
M[n][n] = 1/A[n][n]
for k in range(n-1,0,-1):
M[n][k] = -sum(A[i][k]*M[n][i]
for i in range(n,k,-1))/A[k][k]
return M
Inverse_Bell_Matrix(factorial, 7)
♦ The Bell inverse of a sequence
Given an invertible Bell matrix B with generator g then we call the generator of the inverse of B the Bell inverse of g. Thus we are looking here at a sequence to sequence transformation.
With Sage we can implement this as:
# Given a sequence f returns the inverse Bell sequence of f.
def Inverse_Bell_Sequence(f, dim):
M = Inverse_Bell_Matrix(f, dim)
return [M[n][1] for n in (1..dim-1)]
Some examples:
# [1, -1, 1, -1, 1, -1, 1, -1, ...] maps to # [1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880, ...] print Inverse_Bell_Sequence(lambda n: (-1)^n, 9) # [1, 2, 2, 0, 0, 0, 0, 0, 0, ...] maps to # [1, -2, 10, -80, 880, -12320, 209440, -4188800, ...] print Inverse_Bell_Sequence(lambda n: [1,2,2][n] if n<3 else 0, 9) # [1, 2, 1, 0, 0, 0, 0, 0, 0, ...] maps to # [1, -2, 11, -100, 1270, -20720, 413000, -9726640, ...] print Inverse_Bell_Sequence(lambda n: [1,2,1][n] if n<3 else 0, 9)
More general: Let g be a 0-based sequence then we require g(0) = 0 and g(1) != 0 and write a = g(1), b = g(2), c = g(3), d = g(4) to simplify the notation. The Bell inverse of g then starts
1/a, -b/a^3, (-c*a+3*b^2)/a^5, -(d*a^2-10*b*c*a+15*b^3)/a^7, ...
If we assume a = 1, which will be often the case, this reduces to:
1, -b, -c+3*b^2, -d+10*b*c-15*b^3, ...
| 1 |
| -b |
| 3*b^2 - c |
| 15*b^3 + 10*b*c - d |
| 105*b^4 - 105*b^2*c + 10*c^2 + 15*b*d - e |
Therefore the skeleton of the coefficients is
| 1 |
| -1 |
| 3, -1 |
| -15, 10, -1 |
| 105, -105, 10, 15, -1 |
Looking up this table in the OEIS we find Wolfdieter Lang's A176740 called 'Inversion of e.g.f. formal power series' and his paper E.g.f. Lagrange inversion.
There are two differences: Lang's table misses the leading '1' and the order of coefficients is different. For our purposes the latter is of no importance.
♦ The coefficients of the Bell polynomials
But why stop half way? So let's look at all the coefficients of the inverse Bell polynomials. The coefficients in the table above appear here in column 1.
| n\k | 0 | 1 | 2 | 3 | 4 | 5 |
| 0 | [1] | |||||
| 1 | [0] | [1] | ||||
| 2 | [0] | [-1] | [1] | |||
| 3 | [0] | [3, -1] | [-3] | [1] | ||
| 4 | [0] | [-15, 10, -1] | [15, -4] | [-6] | [1] | |
| 5 | [0] | [105, -105, 10, 15, -1] | [-105, 60, -5] | [45, -10] | [-10] | [1] |
It seems somewhat surprising that this table was not yet in the database. Now it is A268442. Let's also add the table of the coefficients of the Bell polynomials A268441 (also A036040 with different order and column 0 missing).
| n\k | 0 | 1 | 2 | 3 | 4 | 5 |
| 0 | [1] | |||||
| 1 | [0] | [1] | ||||
| 2 | [0] | [1] | [1] | |||
| 3 | [0] | [1] | [3] | [1] | ||
| 4 | [0] | [1] | [3, 4] | [6] | [1] | |
| 5 | [0] | [1] | [10, 5] | [15, 10] | [10] | [1] |
There is a simple check of the tables: replacing the lists by their sums reduces the triangles to the Stirling number tables: to the Stirling numbers of first kind and to the Stirling numbers of the second kind respectively.
The algorithm for building these tables is given below, implemented with Sage.
def bell_polynomial_matrix(dim, inverse):
def bell_polynomial(n):
X = var(['x'+str(i) for i in (0..dim)])
@cached_function
def T(n, k):
if k==0: return k^n
return sum(binomial(n-1,j-1)*T(n-j,k-1)*X[j-1]
for j in (0..n-k+1)).expand()
return [T(n,k) for k in (0..n)]
A = [[f for f in bell_polynomial(n)] for n in range(dim)]
if not inverse: return A
M = [[0 for k in (0..n)] for n in range(dim)]
for n in range(dim):
M[n][n] = 1/A[n][n]
for k in range(n-1,-1,-1):
M[n][k] = expand(-sum(A[i][k]*M[n][i]
for i in range(n,k,-1))/A[k][k])
return M
def coefficient_matrix(M):
def coefficient(p):
c = SR(p).fraction(ZZ).numerator().coefficients()
return [0] if not c else c
return [[coefficient(p) for p in M[n]] for n in range(len(M))]
The two tables above are then generated by the function call
M = bell_polynomial_matrix(8, inverse=true) coefficient_matrix(M)
respectively by the call
L = bell_polynomial_matrix(8, inverse=false) coefficient_matrix(L)
♦ Bell numbers of higher order
Consider the sequence
S0 → T0 → S1 → T1 → S2 → T2 → ...
Here Sn → Tn indicates the Bell transform mapping a sequence Sn to the triangle Tn and Tn → Sn+1 the operator associating a triangle with the sequence of its row sums. For example if we start this iteration with the sequence <1,1,1,...> we get
S0 = A000012 = <1,1,1,...>
T0 = A048993 # Stirling subset numbers,
S1 = A000110 # Bell numbers,
T1 = A264428 # Bell transform of Bell numbers,
S2 = A187761 # second-order Bell numbers,
T2 = A264430 # Bell transform of second-order Bell numbers,
S3 = A264432 # third-order Bell numbers.
Similarly if we start the iteration with the sequence <1,-1,1,-1,...> we get the higher order complementary Bell numbers.
| Sequence | Bell transform |
Inverse Bell transform |
|
| [1,1,1,1,...] | A000012 | A048993 | A048994 |
| Bell numbers | A000110 | A264428 | A264429 |
| Second order Bell numbers | A187761 | A264430 | A264431 |
| Third order Bell numbers | A264432 | A264433 | A264434 |
| [1,-1,1,-1,...] | A033999 | A080417 | A132393 |
| Complementary Bell numbers | A000587 | A264435 | A264436 |
| Second order complementary Bell numbers |
A265023 |
A nice observation is that Joerg Arndt's award winning A187761 (number of maps f: [n] → [n] with f(x) ≤ x and f(f(x)) = f(f(f(x)))) here turns up as the second-order Bell numbers.
def bell_second_order(generator, n):
G = [generator(k) for k in range(n)]
row = lambda n: bell_transform(n, G)
S = [sum(row(k)) for k in range(n)]
return bell_transform(n, S)
Arndt's sequence can now be easily computed from the constant sequence 1,1,1,... as:
[sum(s) for s in [bell_second_order(lambda k: 1, n) for n in range(10)]]
The sequence is also column 1 in the triangle A264430.
[bell_second_order(bell_number, n) for n in range(8)] [[1], [0, 1], [0, 1, 1], [0, 2, 3, 1], [0, 6, 11, 6, 1], [0, 23, 50, 35, 10, 1], [0, 106, 268, 225, 85, 15, 1], [0, 568, 1645, 1603, 735, 175, 21, 1]]
♦ Associated Stirling subset numbers
As we have seen in the table above the Stirling subset numbers are generated by applying the Bell transform to the simplest sequence of positive numbers: the all 1's sequence.
On the other hand one of the simplest methods to transform a sequence is the Hilbert's Hotel Transform (HHT) also known as the 'shift right operation' in CS communities: move simultaneously cell n to cell n+1 and 0 to cell 0.
So what happens if we apply this idea to 1,1,1,... ? We define in Sage:
hht = lambda k: lambda n: 1 if n>=k else 0 for n in (0..5): print bell_matrix(hht(n), 12)
In terms of A-numbers we get for n = 0, 1, 2,... the triangles:
A048993, A008299, A059022, A059023, A059024, A059025.
which are called the (n+1)-associated Stirling subset numbers following Comtet.
However compared to our construction in the last section the construction of this sequence of sequences has a disadvantage: apart from the case n = 0 the resulting matrices are singular; in other words, the inverse Bell transform cannot be applied to these generators.
♦ The Bell transform of multifactorials
An overview about multifactorials can be found in my July 2011 blog post. In Sage parlance the general definition is:
multifactorial = lambda a,b: lambda n: prod(a*k + b for k in (0..n-1))
Here the multifactorials depends on two parameter, a and b. For example multifactorial(2, 2) are the double factorials of even numbers (2*n)!! = 2^n*n!, A000165. Note that multifactorial(2, 2) is a function, evaluated at n=4 gives multifactorial(2, 2)(4) = 384.
Now we can study the Bell transform of multifactorials.
for a in (1..4):
for b in (1..a):
print "Bell matrix generated by multifactorial", (a,b)
bell_matrix(multifactorial(a,b), 6)
print "Inverse Bell matrix generated by multifactorial", (a,b)
inverse_bell_matrix(multifactorial(a,b), 6)
The table below shows the output in terms of A-numbers. From the 20 basic cases 17 have been in the OEIS (the missing three were added by the author). Looking closer at the entries it becomes obvious that the Bell transform of multifactorials was never studied systematically before.
| Type | Multi factorial |
Bell transform |
Inverse Bell transform |
| MF1,1 | A000142 | A132393 | A048993 |
| MF2,1 | A001147 | A132062 | A122848 |
| MF2,2 | A000165 | A039683 | A075497 |
| MF3,1 | A032031 | A203412 | A265605 |
| MF3,2 | A007559 | A004747 | A049404 |
| MF3,3 | A008544 | A051141 | A075498 |
| MF4,1 | A007696 | A265606 | A265604 |
| MF4,2 | A001813 | A119274 | A119275 |
| MF4,3 | A008545 | A265606 | A049410 |
| MF4,4 | A047053 | A051142 | A075499 |
Some special cases of the Bell-inverse of multi-factorials are:
bell_inverse(MFn, n) = n^k for k≥0.
bell_inverse(MFn, n-1) = falling_factorial(n - 1, k) for k≥0.
bell_inverse(MFn, 1) = MFn-1, n-2 (for n≥3).
bell_inverse(MF2*n, n) = [1, n, 0, 0, 0,...].
bell_inverse(MF2*n, 2*n-2) = A161381_row(n-1).
♦ The Bell transform of powers
The power functions can be seen as the degenerated cases of multifactorials: as those cases where a = 0 in the definition of multifactorials.
multifactorial = lambda a,b: lambda n: prod(a*k + b for k in (0..n-1))
These cases have been studied by Wolfdieter Lang. Lang calls the Bell transform of powers "Stirling2 triangles with scaled diagonals".
for n in range(6):
bell_matrix(lambda k: n^k, 7)
Lang calls the signed inverse Bell transform of powers "Generalized Stirling number triangles of first kind".
for n in (0..5):
bell_matrix(lambda k: n^k, 7).inverse()
| n | Powers of n |
Bell transform |
Inverse Bell transform |
| 0 | A000007 | A023531 | A023531 |
| 1 | A000012 | A048993 | A132393 |
| 2 | A000079 | A075497 | A039683 |
| 3 | A000244 | A075498 | A051141 |
| 4 | A000302 | A075499 | A051142 |
| 5 | A000351 | A075500 | A051150 |
| 6 | A000400 | A075501 | A051151 |
| 7 | A000420 | A075502 | A051186 |
| 8 | A001018 | A075503 | A051187 |
| 9 | A001019 | A075504 | A051231 |
The Bell-inverse of the power functions are
bell_inverse(n^k) = n^k*k! for k≥0.
♦ The Bell transform of rising factorials
risingfactorial = lambda n: lambda k: rising_factorial(n, k)
for n in (0..7):
print bell_matrix(risingfactorial(n), 7)
print inverse_bell_matrix(risingfactorial(n), 7)
The Bell transform of the rising factorial and its inverse lead to:
| n | Rising factorial |
Bell transform |
Inverse Bell transform |
| 0 | A000007 | A023531 | A023531 |
| 1 | A000142 | A132393 | A048993 |
| 2 | A000142* | A111596 | A111596 |
| 3 | A001710 | A046089 | A035342 |
| 4 | A001715 | A049352 | A035469 |
| 5 | A001720 | A049353 | A049029 |
| 6 | A001725 | A049374 | A049385 |
| 7 | A001730 | A134141 | A092082 |
The Bell-inverse of rising factorials are multifactorials:
bell_inverse(risingfactorial(n)) = MF(n-1,1)(k+1) for k≥0.
♦ The Bell transform of falling factorials
Similarly for the Bell transform of falling factorials:
fallingfactorial = lambda n: lambda k: falling_factorial(n, k)
for n in (0..5):
print bell_matrix(fallingfactorial(n), 7)
print inverse_bell_matrix(fallingfactorial(n), 7)
| n | Falling factorial A008279_row |
Bell transform |
Inverse Bell transform |
| 0 | row(0) | A023531 | A023531 |
| 1 | row(1) | A049403 A122848 |
A132062 |
| 2 | row(2) | A049404 | A004747 |
| 3 | row(3) | A049410 | A000369 |
| 4 | row(4) | A049424 | A011801 |
| 5 | row(5) | A049411 | A013988 |
The Bell-inverse of falling factorials are multifactorials:
bell_inverse(fallingfactorial(n)) = MF(n+1,n)(k) for k≥0.
♦ The Bell transform of monotonic factorials
There seems to be no standard name for the function
monotonicfactorial = lambda r: lambda n: rising_factorial(r, n)/factorial(n)
I therefore dubbed it monotonic factorial for the purpose of this post (as it counts the number of monotone words with length n over an r-alphabet; most often it is written as binomial(n+r-1, n) ).
for n in (0..4):
print bell_matrix(monotonicfactorial(n), 7)
print inverse_bell_matrix(monotonicfactorial(n), 7)
Astonishing little is found in the OEIS about this class of sequences, with the notable exception of Roger L. Bagula's A137452 (Olivier Gérard's A061356) which is the inverse Bell transform of A000169. Certainly more triangles of this type should be entered into the OEIS.
♦ An implementation for Maple
The Maple implementation of the Bell matrix given below is compatible with the interface and conventions of Sloane's Maple transforms.
# BELL: Computes the Bell matrix of a sequence.
# Given the list [a(0),...,a(dim)] returns the Bell matrix
# with dimension dim+1 of the sequence a.
BELL := proc(a) local M, dim, n, k:
if whattype(a) <> list then RETURN([]) fi:
dim := nops(a);
M := Matrix(dim, shape=triangular[lower]); M[1,1] := 1;
for n from 1 to dim-1 do M[n+1,2] := a[n] od;
for n from 1 to dim-1 do for k from 1 to n do
M[n+1,k+1] := add(binomial(n-1,j-1)*M[n-j+1,k]*M[j+1,2],
j=1..n-k+1) od od;
RETURN(M) end:
BELLi:= proc(a)
if whattype(a) <> list then RETURN([]) fi:
RETURN(linalg[inverse](BELL(a))) end:
# Examples:
a := [seq(1,n=0..8)]: BELL(a); # A048993
b := [seq(n!,n=0..8)]: BELL(b); # A132393
c := [seq(n,n=1..9)]: BELL(c); # A059297
d := [seq((-1)^n*(n+1)^n,n=0..8)]: BELL(d); # A137452
e := [seq(`if`(n::odd,0,2*n!),n=0..8)]: BELL(e); # A137513
f := [seq(`if`(n<2,(-1)^n,0),n=0..8)]: BELL(f); # A104556
g := [seq(doublefactorial(2*n-1),n=0..8)]: BELL(g); # A001497
♦ Index of triangles generated by the Bell transform
The format of the list is [An[Ak]] which means that the triangle An is generated by the sequence Ak. If Ak is not in the OEIS than Ak might point to a related sequence or to dummy sequence Axxxxxx.
[A000369[A008545]] [A001497[A001147]] [A004747[A008544]] [A008275[A133942]] [A008277[A000012]] [A008296[A265313]] [A008297[A133942]] [A008298[A038048]] [A011801[A008546]] [A013988[A008543]] [A023531[A000007]] [A035342[A001147]] [A035469[A007559]] [A038455[A006963]] [A039621[A177885]] [A039683[A000165]] [A039692[A039647]] [A039810[A000110]] [A039811[A000258]] [A039812[A000307]] [A039813[A000357]] [A039814[A003713]] [A039815[A000268]] [A039816[A000310]] [A039817[A000359]] [A046089[A001710]] [A048176[A051262]] [A048786[A034177]] [A048993[A000012]] [A048994[A133942]] [A049029[A007696]] [A049218[A005359]] [A049352[A001715]] [A049353[A001720]] [A049374[A001725]] [A049385[A008548]] [A049403[A000161]] [A049404[A000925]] [A049410[A004552]] [A049411[A008279]] [A049424[A265609]] [A051141[A032031]] [A051142[A047053]] [A051150[A052562]] [A051151[A047058]] [A051186[A051188]] [A051187[A051189]] [A051231[A051232]] [A059297[A000027]] [A059298[A000027]] [A059419[A009006]] [A060281[A001865]] [A061356[A000169]] [A075497[A000079]] [A075498[A000244]] [A075499[A000302]] [A075500[A000351]] [A075501[A000400]] [A075502[A000420]] [A075503[A001018]] [A075504[A001019]] [A075505[A011557]] [A075525[A265024]] [A078521[A038048]] [A079621[A002866]] [A079638[A029767]] [A079639[A006252]] [A079640[A007840]] [A079641[A000629]] [A079642[A089064]] [A086915[A052849]] [A088729[A000670]] [A088814[A000262]] [A092082[A008542]] [A104556[A003475]] [A105278[A000142]] [A105599[A000272]] [A105786[A000272]] [A105819[A055860]] [A106239[A057500]] [A111246[A130716]] [A111593[A155585]] [A111594[A005359]] [A111596[A155456]] [A119274[A001813]] [A119275[A130706]] [A121408[A177145]] [A122848[Axxxxxx]] [A122850[A001147]] [A125553[A208529]] [A129062[A000629]] [A130123[A000038]] [A130191[A000110]] [A130534[A000142]] [A131222[A029767]] [A132056[A045754]] [A132062[A001147]] [A132393[A000142]] [A134141[A001730]] [A135338[A155456]] [A135494[A153881]] [A136590[A136591]] [A136595[A048287]] [A136630[A000035]] [A136656[A155456]] [A137312[A052849]] [A137320[A066459]] [A137339[A052560]] [A137378[A052510]] [A137431[Axxxxxx]] [A137433[Axxxxxx]] [A137452[A213236]] [A137513[Axxxxxx]] [A143395[A000225]] [A143543[A001187]] [A144402[A000161]] [A144633[A144636]] [A144634[A144636]] [A144644[A001899]] [A145520[A000040]] [A147308[A027641]] [A147309[A022902]] [A147311[A122045]] [A147312[A122045]] [A147315[A000111]] [A151509[Axxxxxx]] [A151511[Axxxxxx]] [A166317[A002436]] [A166318[A002436]] [A171996[Axxxxxx]] [A171998[Axxxxxx]] [A174893[A000142]] [A184962[A000670]] [A185285[A004123]] [A185296[Axxxxxx]] [A185415[Axxxxxx]] [A185419[A143523]] [A185422[A080635]] [A185690[A056594]] [A185951[A193356]] [A186366[A000111]] [A187082[Axxxxxx]] [A187084[Axxxxxx]] [A188062[Axxxxxx]] [A188066[Axxxxxx]] [A188832[A009843]] [A189898[A003027]] [A191249[A062738]] [A194938[A039647]] [A195204[A076726]] [A195205[Axxxxxx]] [A202183[Axxxxxx]] [A202184[Axxxxxx]] [A202185[Axxxxxx]] [A202189[Axxxxxx]] [A202190[Axxxxxx]] [A203412[A007559]] [A209849[A155585]] [A215771[A001710]] [A215861[A215851]] [A217756[A129271]] [A223511[A045755]] [A223512[A045756]] [A223513[A045757]] [A223514[Axxxxxx]] [A223515[Axxxxxx]] [A223516[Axxxxxx]] [A223517[Axxxxxx]] [A223518[Axxxxxx]] [A223522[Axxxxxx]] [A225171[A225170]] [A227450[A007395]] [A228534[A009444]] [A228550[A033678]] [A228859[A001832]] [A228892[A002027]] [A247232[A001929]] [A256041[A005212]] [A256042[A256033]] [A256892[A000262]] [A256893[A000670]] [A259286[Axxxxxx]] [A264428[A000110]] [A264429[Axxxxxx]] [A264430[A125273]] [A264431[Axxxxxx]] [A264433[A047889]] [A264434[A011634]] [A264435[A000587]] [A264436[A007549]] [A265314[Axxxxxx]] [A265602[Axxxxxx]] [A265608[A000262]]
In particular this table shows that the coefficients of the Touchard or Bell polynomials (A048993), the Abel polynomials (A137452), the Mittag-Leffler polynomials (A137513), the modified Hermite polynomials (A104556) and the Bessel polynomials (A001497) all are Bell matrices generated by elementary sequences.
♦ References
- D. E. Knuth, Convolution polynomials, Mathematica J. 2.1 (1992), no. 4, 67-78.
Accompanying, there is a SageMath Jupyter Notebook. It can be read and downloaded from GitHub.
