This site is supported by donations to The OEIS Foundation.

User:Peter Luschny/AignerTriangles

From OeisWiki

Jump to: navigation, search

Aigner Triangles

KEYWORDS: Orthogonal polynomials, Jacobi fractions, generalized Motzkin numbers, weighted Motzkin paths, binary sequence operation.

Concerned with sequences:

A000165, A000296, A000698, A001006, A023443, A026300, A049310, A049347, A053121, A056594, A064189, A066325, A098742, A099174, A104562, A123023, A123023, A126120, A137286, A137338, A216916, A217537, A217538.

Contents

An operation ZN x ZN→ ZN

Let us denote the nonnegative numbers by N and the integers by Z. Given two sequences

  • s: N → Z
  • t: N → Z

we associate a lower triangular matrix A, build recursively by

  • A(0,0) = 1; A(0,k) = 0;
  • A(n,k) = A(n-1,k-1)+s(k)*A(n-1,k)+t(k+1)*A(n-1,k+1).

This matrix will be called the Aigner triangle of s and t provided t(0) = 1 and t(n) <> 0 for all n.

Let dim denote the number of rows of the triangle. Than we can compute the Aigner triangle of s and t with Sage as

def ATri(s,t,dim):
    T = matrix(SR,dim,dim)
    for n in range(dim): T[n,n] = 1
    for n in (1..dim-1):
        for k in (0..n-1):
            T[n,k] = T[n-1,k-1]+s(k)*T[n-1,k]+t(k+1)*T[n-1,k+1]
    return T

Note that we have by definition always the diagonal equal to 1. What we are most interested in is the first column of this triangle.

def ASeq(s,t,dim): return ATri(s,t,dim).column(0)

Thus we have an operation ZN x ZN → ZN and we will call the resulting sequence the generalized Motzkin numbers corresponding to s, t.

Modified Charlier polynomials

... and indecomposable set partitions without singletons

A first example of this operation: Let

s = lambda n: n+1
t = lambda n: n+1

then ATri(s, t, 9) gives

    1     
    1     1     
    3     3     1     
    9    12     6     1     
   33    51    34    10     1
  135   237   193    79    15     1
  609  1188  1132   584   160    21     1
 2985  6381  6920  4268  1510   293    28     1
15747 36507 44213 31542 13576  3464   497    36     1

and ASeq(s, t, 12) gives the generalized Motzkin numbers corresponding to (s, t) = (n+1, n+1), which are

1, 1, 3, 9, 33, 135, 609, 2985, 15747, 88761, 531561, ...

Looking things up in the database we find the number of indecomposable set partitions of [1..n+2] without singletons (in our enumeration), Knuth's A098742, and the triangle A216916.

Since the Aigner triangle of s and t is invertible we may also consider the inverse matrix, which is again a lower triangular array with diagonal equal to 1.

def IATri(s,t,dim):
    M = ATri(s,t,dim)
    return M.inverse()

For example IATri(s, t, 9) with the above s and t gives

     1      
    -1      1
     0     -3      1
     3      6     -6      1
   -12     -9     26    -10      1
    45      3   -109     71    -15      1
  -198     81    501   -475    155    -21     1
  1071   -786  -2663   3329  -1455    295   -28    1
 -6984   6711  16510 -25495  13729  -3647   511  -36   1

This triangle has a nice interpretation: it gives the coefficients of the modified Charlier polynomials  ordered by rising powers. In the database they were entered by R. L. Bagula, A137338.

By the way, the modified Charlier polynomials satisfy the recurrence

C(n, x) = if n > 0 then (x-n)*C(n-1, x) - n*C(n-2, x)
          elif n = 0 then 1 else 0. 

We call the first column of the inverse of an Aigner triangle the inverse generalized Motzkin numbers corresponding to (s, t). They can be computed by

def IASeq(s,t,dim): return IATri(s,t,dim).column(0)

In our example we find:

1, -1, 0, 3, -12, 45, -198, 1071, -6984, 53217, -462330, ..

This sequence is not yet in the database.

Orthogonal polynomials

The way we proceeded to compute the inverse triangle was laborious. First we computed the triangle ATri(s, t) and then computed the inverse. However, there is a much better way to compute IATri(s, t), namely again by recursion directly from s and t.

  • A*(0,0) = 1; A*(0,k) = 0;
  • A*(n+1,k) = A*(n,k-1)-s(n)A*(n,k)-t(n)A*(n-1,k).

Now we can replace the routine IATri by the more efficient one:

def IATri(s,t,dim):
    T = matrix(SR,dim,dim)
    for n in range(dim): T[n,n] = 1
    for n in (0..dim-2):
        for k in (0..n):
            T[n+1,k] = T[n,k-1]-s(n)*T[n,k]-t(n)*T[n-1,k]
    return T

We have already seen above that we could associate a family of polynomials (the modified Charlier polynomials) with the Aigner triangle corresponding to (s(n)=n+1, t(n)=n+1). Such an association is possible in general. To this end we define a recursive system

  • p0(x) = 1;
  • pn+1(x) = (x-sn)pn(x) - tnpn-1(x).                  (*)

Implemented with Sage:

def Poly(s, t, n, x):
    if n < 0 : return 0
    if n == 0 : return 1
    p = (x-s(n-1))*Poly(s,t,n-1,x)-t(n-1)*Poly(s,t,n-2,x)
    return expand(p)

Continuing our example:

for n in (0..7): Poly(s,t,n,x)

1
x   -  1
x^2 -  3*x
x^3 -  6*x^2 +   6*x   +    3
x^4 - 10*x^3 +  26*x^2 -    9*x   -   12
x^5 - 15*x^4 +  71*x^3 -  109*x^2 +    3*x   + 45
x^6 - 21*x^5 + 155*x^4 -  475*x^3 +  501*x^2 + 81*x - 198

And these are again the modified Charlier polynomials, as we expected.

Now is there something special about polynomials as defined in (*) above? Indeed there is. It is described in a theorem due to Favard:

A sequence of polynomials pn(x) forms an orthogonal system if and only if (*) holds for pn(x) for some sequences s and t.

The associated Hankel matrix

To prove this theorem we have to introduce two more entities related with every Aigner triangle of s and t.

The first is trivial: We associate with t the sequence of partial products Tn = t0t1t1...tn and for notational convenience the diagonal matrix T, hinted by:

T_0               
      T_1           0
            T_2
 0            .
                    .

T can be easily computed with Sage:

def T(t, dim):
    M = matrix(SR, dim, dim)
    for n in range(dim): 
        M[n,n] = mul(t(i) for i in (0..n))
    return M

Continuing our example we find:

T(t,7)    

[   1    0    0    0    0    0    0]
[   0    2    0    0    0    0    0]
[   0    0    6    0    0    0    0]
[   0    0    0   24    0    0    0]
[   0    0    0    0  120    0    0]
[   0    0    0    0    0  720    0]
[   0    0    0    0    0    0 5040]

The second entity which we want to associate with an Aigner triangle is much more gripping: it is the Hankel matrix of ASeq(s, t), the generalized Motzkin numbers corresponding to s, t (see above).

def Hankel(f, dim):
    H = matrix(SR, dim, dim)
    for i in range(dim):
        for j in range(dim):
            H[i,j] = f[i+j]
    return H      
    
def Han(s, t, dim):
    b = ASeq(s,t,2*dim-1) # generalized Motzkin numbers wrt s, t.
    return Hankel(b, dim)  

Now the technical bottom line can be simply stated as:

det Hn = det Tn for all n.

Let's check this with our example:

[Han(s,t,n).det() for n in (1..7)] 
[1, 2, 12, 288, 34560, 24883200, 125411328000]

[T(t,n).det() for n in (1..7)] 
[1, 2, 12, 288, 34560, 24883200, 125411328000]

These are the superfactorials A000178.

The general bottom line can now be stated as:

A sequence of numbers are generalized Motzkin numbers corresponding to some s, t if and only if det Hn <> 0 for all n.

We will also look at the row sums and the alternating row sums of the triangles:

def Sum(L): 
    return add(L[k] for k in range(len(L)))
def AltSum(L): 
    return abs(add((-1)^k*abs(L[k]) for k in range(len(L))))

Let us summarize our findings so far:

s(n) = n+1,  t(n) = n+1
  Triangle Sequence Sum AlterSum
AT A216916 A098742 - -
AT-1 A137338 - - -
  Hankel: A000178 Poly: modified Charlier


Motzkin numbers and shifted Chebyshev polynomials

Now let us look at a second example, which is the most characteristic one: setting s(n) = t(n) = 1. This will give us the Urform of all the relations studied here.

s = lambda n: 1
t = lambda n: 1

ATri(s, t, 9)     # A064189 A026300
[  1   0   0   0   0   0   0   0   0]
[  1   1   0   0   0   0   0   0   0]
[  2   2   1   0   0   0   0   0   0]
[  4   5   3   1   0   0   0   0   0]
[  9  12   9   4   1   0   0   0   0]
[ 21  30  25  14   5   1   0   0   0]
[ 51  76  69  44  20   6   1   0   0]
[127 196 189 133  70  27   7   1   0]
[323 512 518 392 230 104  35   8   1]

ASeq(s,t,12)      # A001006 
1, 1, 2, 4, 9, 21, 51, 127, 323, 835, 2188, 5798

ASum(s,t,12)      # A005773
1, 2, 5, 13, 35, 96, 267, 750, 2123, 6046, 17303

IATri(s, t, 9)  # A104562
[  1   0   0   0   0   0   0   0   0]
[ -1   1   0   0   0   0   0   0   0]
[  0  -2   1   0   0   0   0   0   0]
[  1   1  -3   1   0   0   0   0   0]
[ -1   2   3  -4   1   0   0   0   0]
[  0  -4   2   6  -5   1   0   0   0]
[  1   2  -9   0  10  -6   1   0   0]
[ -1   3   9 -15  -5  15  -7   1   0]
[  0  -6   3  24 -20 -14  21  -8   1]

IASeq(s,t,12)  # A049347
1, -1, 0, 1, -1, 0, 1, -1, 0, 1, -1, 0

IASum(s,t,12)  # A056594
1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0, 

[Han(s,t,n).det() for n in (1..7)] 
[1, 1, 1, 1, 1, 1, 1]   # A000012 

So we can summarize: The sequence pair s(n) = t(n) = 1 maps to the Motzkin triangle and the Motzkin numbers.  The row sums are the number of directed animals of size n A005773 and the alternating sums the Riordan numbers A005043.

s(n) = 1 ,  t(n) = 1
  Triangle Sequence Sum AlterSum
AT A064189, A026300 A001006 A005773 A005043
AT-1 A104562 A049347 A056594 -
  Hankel: A000012 Poly: Chebyshev


Catalan numbers and Chebyshev polynomials

Our third example gives the most important special case of the Aigner triangles. Here we set s(n) = 0 and t(n) = 1.

s = lambda n: 0
t = lambda n: 1

ATri(s, t, 9)  # A053121
[ 1  0  0  0  0  0  0  0  0]
[ 0  1  0  0  0  0  0  0  0]
[ 1  0  1  0  0  0  0  0  0]
[ 0  2  0  1  0  0  0  0  0]
[ 2  0  3  0  1  0  0  0  0]
[ 0  5  0  4  0  1  0  0  0]
[ 5  0  9  0  5  0  1  0  0]
[ 0 14  0 14  0  6  0  1  0]
[14  0 28  0 20  0  7  0  1]

ASeq(s,t,12)  # A126120 
1, 0, 1, 0, 2, 0, 5, 0, 14, 0, 42, 0

ASum(s,t,12)  # A001405
1, 1, 2, 3, 6, 10, 20, 35, 70, 126,

IATri(s, t, 9)  # A049310
[  1   0   0   0   0   0   0   0   0]
[  0   1   0   0   0   0   0   0   0]
[ -1   0   1   0   0   0   0   0   0]
[  0  -2   0   1   0   0   0   0   0]
[  1   0  -3   0   1   0   0   0   0]
[  0   3   0  -4   0   1   0   0   0]
[ -1   0   6   0  -5   0   1   0   0]
[  0  -4   0  10   0  -6   0   1   0]
[  1   0 -10   0  15   0  -7   0   1]

IASeq(s,t,12)   # A056594 
1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0, 

IASum(s,t,12)   # A010892 
1, 1, 0, -1, -1, 0, 1, 1, 0, -1, -1, 0, 

[Han(s,t,n).det() for n in (1..7)] 
[1, 1, 1, 1, 1, 1, 1]  # A000012 

Thus we find: The sequence pair s(n) = 0 and t(n) = 1 maps to the Catalan triangle and the aerated Catalan numbers. The alternating row sums of the inverse triangle are the Fibonacci numbers A000045.

s(n) = 0,  t(n) = 1
  Triangle Sequence Sum AlterSum
AT A053121 A126120 A001405 A001405
AT-1 A049310 A056594 A010892 A000045
  Hankel: A000012 Poly: Chebyshev


Set partitions without singletons

We come back to our first example, however with a small modification.

s = lambda n: n
t = lambda n: n if n > 0 else 1

ATri(s, t, 9)  # A217537
[   1    0    0    0    0    0    0    0    0]
[   0    1    0    0    0    0    0    0    0]
[   1    1    1    0    0    0    0    0    0]
[   1    4    3    1    0    0    0    0    0]
[   4   11   13    6    1    0    0    0    0]
[  11   41   55   35   10    1    0    0    0]
[  41  162  256  200   80   15    1    0    0]
[ 162  715 1274 1176  595  161   21    1    0]
[ 715 3425 6791 7182 4361 1526  294   28    1]

ASeq(s,t,12)  # A000296
1, 0, 1, 1, 4, 11, 41, 162, 715, 3425, 17722, 98253,

ASum(s,t,12)  # A217924 
1, 1, 3, 9, 35, 153, 755, 4105, 24323, 155513, 

IATri(s, t, 9)  
[    1     0     0     0     0     0     0     0     0]
[    0     1     0     0     0     0     0     0     0]
[   -1    -1     1     0     0     0     0     0     0]
[    2    -1    -3     1     0     0     0     0     0]
[   -3     8     5    -6     1     0     0     0     0]
[    4   -31     0    25   -10     1     0     0     0]
[   -5   119   -56   -95    70   -15     1     0     0]
[    6  -533   455   364  -455   154   -21     1     0]
[   -7  2904 -3326 -1428  3059 -1428   294   -28     1]

IASeq(s,t,12)   # (-1)^(n+1)*A023443(n)  
1, 0, -1, 2, -3, 4, -5, 6, -7, 8, -9, 10, -11, 12, -13, 

IASum(s,t,12)   # A165900 
1, 1, -1, -1, 5, -11, 19, -29, 41, -55, 71, -89, 109, 

[Han(s,t,n).det() for n in (1..7)]    # A000178 
[1, 1, 2, 12, 288, 34560, 24883200]  
s(n) = n,  t(n) = n
  Triangle Sequence Sum AlterSum
AT A217537 A000296 A217924 A000012
AT-1 - A023443 A165900 -
  Hankel: A000178 Poly: P


Feynman diagrams and modified Hermite polynomials

Next let's see what happens if we make the choice s(n) = 0 and t(n) = n+1.

s = lambda n: 0
t = lambda n: n+1

ATri(s, t, 9)  # A217538 
[   1    0    0    0    0    0    0    0    0]
[   0    1    0    0    0    0    0    0    0]
[   2    0    1    0    0    0    0    0    0]
[   0    5    0    1    0    0    0    0    0]
[  10    0    9    0    1    0    0    0    0]
[   0   37    0   14    0    1    0    0    0]
[  74    0   93    0   20    0    1    0    0]
[   0  353    0  193    0   27    0    1    0]
[ 706    0 1125    0  355    0   35    0    1]

Aseq(s,t,13)  # aerated A000698
1, 0, 2, 0, 10, 0, 74, 0, 706, 0, 8162, 0, 110410,

Asum(s,t,13)  # 
1, 1, 3, 6, 20, 52, 188, 574, 2222, 7626, 31350,

IATri(s, t, 9)  # A137286
[   1    0    0    0    0    0    0    0    0]
[   0    1    0    0    0    0    0    0    0]
[  -2    0    1    0    0    0    0    0    0]
[   0   -5    0    1    0    0    0    0    0]
[   8    0   -9    0    1    0    0    0    0]
[   0   33    0  -14    0    1    0    0    0]
[ -48    0   87    0  -20    0    1    0    0]
[   0 -279    0  185    0  -27    0    1    0]
[ 384    0 -975    0  345    0  -35    0    1]

IASeq(s,t,12)   # aerated A000165 
1, 0, -2, 0, 8, 0, -48, 0, 384, 0, -3840, 0

IASum(s,t,13)  # 
1, 1, -1, -4, 0, 20, 20, -120, -280, 800, 3600

[Han(s,t,n).det() for n in (1..7)]  # A000178(n+1)
[1, 2, 12, 288, 34560, 24883200, 125411328000]   

The most intriguing observation is the connection with A000698. The inverse triangle is A137286, the coefficients of a variant of the Hermite polynomials. The before missing Aigner triangle is now A217538.

s(n) = 0,  t(n) = n+1
  Triangle Sequence Sum AlterSum
AT A217538 A000698 - -
AT-1 A137286 A000165 - A000932
  Hankel: A000178 Poly: modified-Hermite


Schroeder's third problem and Hermite polynomials

A variant of the foregoing is:

s = lambda n: 0
t = lambda n: n if n > 0 else 1

ATri(s, t, 9)  # A099174 
[  1   0   0   0   0   0   0   0   0]
[  0   1   0   0   0   0   0   0   0]
[  1   0   1   0   0   0   0   0   0]
[  0   3   0   1   0   0   0   0   0]
[  3   0   6   0   1   0   0   0   0]
[  0  15   0  10   0   1   0   0   0]
[ 15   0  45   0  15   0   1   0   0]
[  0 105   0 105   0  21   0   1   0]
[105   0 420   0 210   0  28   0   1]

ASeq(s,t,13)  # A123023 
1, 0, 1, 0, 3, 0, 15, 0, 105, 0, 945, 0, 10395,

ASum(s,t,13)  # A000085 
1, 1, 2, 4, 10, 26, 76, 232, 764, 2620, 9496,

IATri(s, t, 9)  # A066325 
[   1    0    0    0    0    0    0    0    0]
[   0    1    0    0    0    0    0    0    0]
[  -1    0    1    0    0    0    0    0    0]
[   0   -3    0    1    0    0    0    0    0]
[   3    0   -6    0    1    0    0    0    0]
[   0   15    0  -10    0    1    0    0    0]
[ -15    0   45    0  -15    0    1    0    0]
[   0 -105    0  105    0  -21    0    1    0]
[ 105    0 -420    0  210    0  -28    0    1]

IASeq(s,t,12)  # A123023
1, 0, -1, 0, 3, 0, -15, 0, 105, 0, -945, 0, 10395, 

IASum(s,t,13)  # A001464
1, 1, 0, -2, -2, 6, 16, -20, -132, 28, 1216, 936,

[Han(s,t,n).det() for n in (1..7)]  # A000178
[1, 1, 2, 12, 288, 34560, 24883200, 125411328000]   
s(n) = 0, t(n) = n if n > 0 else 1
  Triangle Sequence Sum AlterSum
AT A099174 A123023 A000085 A000085
AT-1 A066325 A123023 A001464 A000085
  Hankel: A000178 Poly: Hermite


Jacobi Fractions

Time to switch back to some notions. The Jacobi fractions are continued fractions of the form


C=\frac{\beta_0 x^2}{1-\alpha_0 x -}\frac{\beta_1 x^2}{1-\alpha_1 x -}\frac{\beta_2 x^2}{1-\alpha_2 x-} \dots

More on the Digital Library of Mathematical Functions: DLMF

Let's raise the question: What happens if we feed our sequence pair s and t into the Jacobi fractions?

# Jcf evaluates a continued fraction of Jacobi type
# using backwards recurrence.
def Jcf(s,t,n):
    c = 0
    for k in range(n,0,-1):
        a = 1-s(k)*x
        b = t(k)*x^2
        c = b/(a-c)
    c = 1/(1-s(0)*x-c)
    return c.simplify_full()
  • The Catalan example (s(n)=0, t(n)=1) gives us:
for n in (0..6): Jcf(s,t,n)
1
-1/(x^2 - 1)
(x^2 - 1)/(2*x^2 - 1)
-(2*x^2 - 1)/(x^4 - 3*x^2 + 1)
(x^4 - 3*x^2 + 1)/(3*x^4 - 4*x^2 + 1)
-(3*x^4 - 4*x^2 + 1)/(x^6 - 6*x^4 + 5*x^2 - 1)
(x^6 - 6*x^4 + 5*x^2 - 1)/(4*x^6 - 10*x^4 + 6*x^2 - 1)

Next let's see what we can recover if we expand these rational functions.

for n in (0..6): n, taylor(Jcf(s,t,n),x,0,2*n)
(0) 1
(1) 1 + 1*x^2  
(2) 1 + 1*x^2 + 2*x^4
(3) 1 + 1*x^2 + 2*x^4 + 5*x^6
(4) 1 + 1*x^2 + 2*x^4 + 5*x^6 + 14*x^8
(5) 1 + 1*x^2 + 2*x^4 + 5*x^6 + 14*x^8 + 42*x^10
(6) 1 + 1*x^2 + 2*x^4 + 5*x^6 + 14*x^8 + 42*x^10 + 132*x^12 
  • In the Motzkin case (s(n)=1, t(n)=1) we find:
for n in (0..5): Jcf(s,t,n)
1/(x - 1)
(x - 1)/(2*x - 1)
-(2*x - 1)/(x^3 + x^2 - 3*x + 1)
-(x^3+x^2-3*x+1)/(x^4-2*x^3-3*x^2+4*x-1)
(x^4-2*x^3-3*x^2+4*x-1)/(4*x^4-2*x^3-6*x^2+5*x-1)
-(4*x^4-2*x^3-6*x^2+5*x-1)/(x^6+2*x^5-9*x^4+10*x^2-6*x+1)

and expanded:
(0) 1
(1) 1+x+2*x^2
(2) 1+x+2*x^2+4*x^3+9*x^4
(3) 1+x+2*x^2+4*x^3+9*x^4+21*x^5+51*x^6
(4) 1+x+2*x^2+4*x^3+9*x^4+21*x^5+51*x^6+127*x^7+323*x^8

Thus we see: Going to infinity with n yields the infinite J-fraction expansion of Sum_{n≥0} ASeqn(s,t) zn. In other words, we found ordinary generating functions of the generalized Motzkin numbers associated with s and t (up to n).

Weighted Motzkin paths

We have seen that a great variety of combinatorial sequences arise from Aigner triangles, which can be computed recursively and have Jacobi fractions as generating functions. They give rise to orthogonal polynomials and indeed (the coefficients of) all possible orthogonal polynomials can be found from Aigner triangles.

But the icing on the cake is that there is a uniform combinatorial interpretation for all these objects. The name is weighted Motzkin paths.

A Motzkin path is a path in R2 from (0, 0) to (n, 0) consisting of a sequence of steps provided that it never falls below the x-axis. Three types of steps are permitted, denoted by (1,1), (1,0) and (1,-1), and called up-steps, horizontal steps and down-steps, respectively.

We will work here with Motzkin paths by letting 1’s represent steps in the direction (1,1) and 0’s represent steps in the direction (1,0) and -1’s represent steps in the direction (1,-1). At least when we are coding. When pretty printing a path we will use instead the letters 'H', 'D' and 'U' for '0', '-1' and '1', respectively.

def steps_to_string(S):
    Z = ""
    for s in S:
        if   s== 0: Z += 'H'  # horizontal step
        elif s==-1: Z += 'D'  # down step
        elif s== 1: Z += 'U'  # up step
    return Z

def PrintAllMotzkinPaths(n):
    P = iter(MotzkinPaths(n))
    for p in P:
        print steps_to_string(p),
    return

PrintAllMotzkinPaths(5)

HHHHH HHHUD HHUHD HHUDH HUHHD HUHDH HUDHH HUDUD 
HUUDD UHHHD UHHDH UHDHH UHDUD UHUDD UDHHH UDHUD
UDUHD UDUDH UUHDD UUDHD UUDDH

Here we used the following elegant Sage function written by Dan Drake to generate the Motzkin paths.

# Copyright (C) 2011 Dan Drake. This program is free software: 
# you can redistribute it and/or modify it under the terms of
# the GNU General Public License as published by the Free 
# Software Foundation, either version 2 of the License, or (at
# your option) any later version. See: http://www.gnu.org/licenses

from sage.combinat.backtrack import GenericBacktracker

class MotzkinPaths(GenericBacktracker):
    def __init__(self, n, h=None):
        GenericBacktracker.__init__(self, [], (0, 0))
        self._n = n
        if h is not None:
            self.max_ht = h
        else:
            self.max_ht = n

    def _rec(self, path, state):

        len, ht = state

        if len < self._n:
            # if length is less than n, we need to keep building 
            # the path, so new length will be 1 longer, and we 
            # don't yield the path yet.
            newlen = len + 1

            # can always tack on an E step, if the new height isn't
            # more than the new number of steps remaining
            if ht <= self._n - newlen:
                yield path + [0], (newlen, ht), False

            # if we're not touching the x-axis, we can yield a path
            # with a downstep at the end
            if ht > 0:
                yield path + [-1], (newlen, ht - 1), False

            # can take an upstep if the new height isn't more than
            # the new number of steps remaining, and if we're not
            # at the maximum height
            if ht + 1 <= self._n - newlen and ht < self.max_ht:
                yield path + [1], (newlen, ht + 1), False
        else:
            # if length is n, set state to None so we stop trying
            # to make new paths, and yield what we've got
            yield path, None, True

# - End of Dan Drake's MotzkinPaths.

Now comes the crucial step. We associate with a path a weight, depending on the function pair s and t of course.

def weight(s,t,P):
    h = 0; w = 1
    for p in P:
        if   p== 1: h += 1
        elif p== 0: w *= s(h)
        elif p==-1: w *= t(h); h -= 1;
    return w

Basically this says: The weight of an up-step is 1, the weight of an horizontal step is given by the s-function and the weight of a down-step by the t-function where the functions are evaluated at the height of the path in this point. The weight of the path is the product of the weight of its steps.

For example we find for Motzkin paths with length 4 and

(I)   s(n) = n+1 and t(n) = n+1 (first example above)
(II)  s(n) = 1 and t(n) = 1     (Motzkin case)
(III) s(n) = 0 and t(n) = 1     (Catalan case)

             [I]  [II]  [III]    
   --------------------------
   w(HHHH) =  1  |  1  |  0
   w(HHUD) =  2  |  1  |  0
   w(HUHD) =  4  |  1  |  0
   w(HUDH) =  2  |  1  |  0
   w(UHHD) =  8  |  1  |  0
   w(UHDH) =  4  |  1  |  0
   w(UDHH) =  2  |  1  |  0
   w(UDUD) =  4  |  1  |  1
   w(UUDD) =  6  |  1  |  1
             --     -     -
             33     9     2     

The bottom line of this construction is:

The sums of the weights over all Motzkin paths of length n are the generalized Motzkin numbers associated with s and t.

Selected bibliography

  • M. Aigner, Catalan and other numbers: a recurrent theme. In: Algebraic Combinatorics and Computer Science, Springer 2001. (What we called Aigner matrix here Aigner calls a 'recursive matrix'. However this term appears with different meanings in the literature.)
  • P. Flajolet, Combinatorial aspects of continued fractions, Discrete Math. 32 (1980), 125-161
  • C. G. J. Jacobi, Werke, vol. VI, p. 385-426.
  • X. Viennot, Une théorie combinatoire des polynômes orthogonaux généraux, Lecture notes, Montréal (1984).
  • H. Wronski, Philosophie de la Technie algorithmique: Loi suprême et universelle. Paris 1815-1817
Personal tools