This site is supported by donations to The OEIS Foundation.

User:Peter Luschny/BellTransform

From OeisWiki

Jump to: navigation, search

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

♦ 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.

Coefficients of inverse Bell polynomials
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).

Coefficients of Bell polynomials
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


Accompanying, there is an Jupyter Notebook which was developed on SageMath Cloud. SageMath is a free open source alternative to Magma, Maple, Mathematica and Matlab. Get a free account on SageMath Cloud, download this notebook, upload it to SMC and start exploring!

Personal tools