The Hitchhiker's Guide to Riordan Squares

Eric Weisstein’s World of Mathematics offers no information on Riordan arrays, the English Wikipedia exactly one sentence. Since many sequences in the OEIS are Riordan arrays and identifying them illuminates their inner structure this is unsatisfactory.

A Riordan array is an infinite lower triangular matrix depending on two generating functions. Riordan squares are a subspecies of the Riordan arrays. Like a Riordan array a Riordan square is an infinite lower triangular matrix, but only depends on one generating function. A Riordan square can often be used as a replacement for a Riordan array which simplify things considerably. To our knowledge, Riordan squares as introduced here are not described elsewhere.

Our approach is informal and leaves out the purely mathematical side as far as possible. Instead, we will look at Riordan squares from the perspective of computation. A familiarity with reading algorithmic presentations will be useful.

On our walk through the examples given in the OEIS we will touch on many other subject areas, in particular continued fractions of Jacobi and Stieltjes type and Deléham's Δ operation. As a fascinating application the construction of orthogonal polynomials using Riordan squares is sketched.

The closing part of the article contains a list of important examples from the field of enumerative combinatorics and links them to the OEIS, to which we added several new entries during the writing of this article.

The transform RS: ZN ↦ Mat[Z]

The definition of a Riordan square is based on a single operation, the Cauchy product (or convolution). Given two lists a = [a0, ..., an] and b = [b0, ..., bn], and denoting the Cauchy product with #, we set

a # b = Sum_{k=0..n} a(k) * b(n - k).

To compute the Riordan square of a sequence S we first set column 0 of the triangle R to the values of S, R(n, 0) = S(n) for n≥0.

Thereafter, i.e. for n≥1, R(n, k) is set to the convolution of the slice [R(1, 0), ..., R(k+1, 0)] of column 0 with the slice [R(k-1, k-1), ..., R(2k-1, k-1)] of column k - 1.

(0, 0)  
(1, 0) ...  
(2, 0) ... ...  
... ... ... (k−1, k−1)  
(k, 0) ... ... (k, k−1) ...  
 (k+1, 0)  ... ... ... ... ...
... ... ... (2k−2, k−1) ... ...
... ... ... (2k−1, k−1) ... ...
... ... ... ...  (n,  k)  ...

R(n, k) = [R(1, 0), ..., R(k+1, 0)] # [R(k-1, k-1), ..., R(2k-1, k-1)]

This formula translates to the function below which computes the Riordan square of S as a lower triangular matrix. Note that indexing starts with n = 0 and k = 0, i.e. the function returns a (0, 0)-based triangle.

    def RiordanSquare(S):
        N = len(S)
        M = matrix(QQ, N, N)
        for n in (0..N-1): M[n, 0] = S[n]
        for k in (1..N-1):
            for m in (k..N-1):
                M[m, k] = sum(M[j, k-1]*S[m-j] for j in (k-1..m-1))
        return M

Assuming the sequence is S = a, b, c, d, e, ... the skeleton of the Riordan square of S starts:

    a;
    b, ab;
    c, ac+b^2,     ab^2;
    d, ad+2bc,     2abc+b^3,            ab^3;
    e, ae+2bd+c^2, 2abd+ac^2+3cb^2,     3ab^2c+b^4,         ab^4;
    f, af+2be+2cd,2abe+2dac+3db^2+3bc^2,3ab^2d+3abc^2+4cb^3,4ab^3c+b^5,ab^5;

As a matrix the triangle is invertible if and only if a and b are not zero. This condition is natural for Riordan squares and assumed in applications of the transform although not needed for the algorithm to work. Particular nice cases arise when a, b ∈ {-1, 1}. In the case a = b = 1 the ordinary RS reduces to:

    1;
    1, 1;
    c, 1+c,      1;
    d, d+2c,     2c+1,           1;
    e, e+2d+c^2, 2d+3c+c^2,      3c+1,       1;
    f, f+2e+2cd, 2e+3d+2cd+3c^2, 3d+3c^2+4c, 4c+1, 1;

The mother of all Riordan squares

If one sets up a new transformation one likes to know on what value the unity is mapped. In our case the unity is the sequence 1, 1, 1, ..., so let’s start with that. Let a = b = c = d = e = f = 1. Then the triangle reduces to

    1;
    1, 1;
    1, 2,  1;
    1, 3,  3,  1;
    1, 4,  6,  4, 1;
    1, 5, 10, 10, 5, 1;

This shows that R(1, 1, 1, ...) is Pascal's triangle; this is good news. We can say that the ordinary Riordan square is a refinement of Pascal's triangle, to which it reduces in the most simple case .

Let's examine the case in the figure above. The left (yellow) stripe consists only of '1's which appear as factors in the Cauchy product. So we do not have to pay attention to them. The values in the green stripe remain, the sum of which gives the value in the blue cell. This is of course one of the well known methods to calculate Pascal's triangle.

Alternating row sums

Let R be a Riordan square matrix with R(0, 0) = 1 then the alternating row sums vanish for all n ≥ 1. There is also some kind of dual: if S(0) = R(0, 0) = -1 then the row sums vanish for n ≥ 1. Following Don Knuth we use here the convention 0^0 = 1.

     0^n = Sum_{k=0..n} (-1)^k*R(n, k)  if S(0) =  1
    -0^n = Sum_{k=0..n} R(n, k)         if S(0) = -1

This also means that if a sequence S with S(0) = 1 generates the triangle R, then S can be recovered from R without the first column through

    S(n) = -Sum_{k=1..n} (-1)^k*R(n, k) for n ≥ 0.

An adapter for series and functions

To start the transform from sequences like we did above is not typical. For historical and mathematical reasons it is more common to start from a series. From our (software-oriented) viewpoint this is just a different interface of the transform for which we have to provide an adapter.

The implementation of this adapter is in the appendix. Applied to the function RiordanSquare the input becomes polymorphic. It does not matter if the input is a finite sequence (list), a generating function or an ordinary function — it will behave the same in all cases.

Just enter what you have at hand. The syntax is in all cases (almost) the same. First give the representation of the sequence then the desired length of the output. If the input is a list the latter is omitted and set to the length of the list. For an example consider the positive integers.

    RiordanSquare([1, 2, 3, 4, 5, 6, 7, 8])
    RiordanSquare(lambda n: n + 1, 8)
    RiordanSquare(1/(1-x)^2, 8)

    1
    2    2
    3    7    4
    4   16   20    8
    5   30   61   52   16
    6   50  146  198  128   32
    7   77  301  575  584  304   64
    8  112  560 1408 1992 1616  704  128

Recall that indexing is 0-based in our setup. Thus if we want to start the first column with 1 we have to give the function n ↦ n + 1 as a parameter.

The exponential case

There is a second flavor of the Riordan squares, the exponential Riordan squares. We implement the exponential case as a post-processor to the ordinary case in a function called ExponentialWeights. This name hints that there are other post-processors or weights which could be applied to the Riordan square, but we leave this topic for another time.

To include the exponential case in our setup we expand the function like this:

    def RiordanSquare(S, expo=false):
        N = len(S)
        M = matrix(QQ, N, N)
        for n in (0..N-1): M[n, 0] = S[n]
        for k in (1..N-1):
            for m in (k..N-1):
                M[m, k] = sum(M[j, k-1]*S[m-j] for j in (k-1..m-1))
        if expo: return ExponentialWeights(M)
        return M

    def ExponentialWeights(M):
        N = M.nrows()
        u = 1
        for k in (1..N-1):
            u *= k
            for m in (0..k):
                j = u if m == 0 else j/m
                M[k, m] *= j
        return M

If you want the exponential version of the algorithm to be applied call the function in the form RiordanSquare(S, expo=true). By default 'expo' is set to 'false' and does not need to be indicated at all.

The skeleton of the exponential Riordan Square triangle Re(S) starts:

    a;
    b,   ab;
    2c,  2ac+2b^2,        ab^2;
    6d,  6ad+12bc,        6abc+3b^3,           ab^3;
    24e, 24ae+48bd+24c^2, 24abd+12ac^2+36cb^2, 12ab^2c+4b^4, ab^4;

The call RiordanSquare([1, 1, 1, 1, 1, 1, 1], expo=true) returns the first seven rows of the coefficients of the unsigned Laguerre polynomials, n! L(n, x).

    1;
    1,   1;
    2,   4,   1;
    6,   18,  9,   1;
    24,  96,  72,  16,  1;
    120, 600, 600, 200, 25, 1;

Here the sequence 1, 1, 1, ... maps to an important family of orthogonal polynomials. The Laguerre polynomials have countless applications in physics and other sciences. So we have also a promising start for studying the exponential case.

In the frequent case a = b = 1 the inverse matrix of the exponential RS starts:

     1;
    -1,  1;
     2, -2-2c,                               1;
    -6, -6d+6c+12c^2+6,                     -3-6c,              1;
    24, -24e+24d-48c^2+120cd-24c-120c^3-24, -24d+24c+60c^2+12, -4-12c, 1;

These triangles play an important role for generating orthogonal polynomials. Their trademark are the alternating factorials in the first column. Examples will be given below.

The product Z[[x]] X Z[[x]] ↦ Mat[Z]

The Riordan product is a map a, b ↦ [a, b] associating two formal power series a, b with a lower triangular matrix [a, b]. The Riordan square is the case a = b of the Riordan product. Formally we can describe the Riordan square as a transform RS: Z[[x]] ↦ Mat[Z] which maps power series over the integers to (lower triangular) integer matrices.

A Riordan array is a lower triangular matrix M which is the Riordan product of two power series a, b where one requires that a(0) ≠ 0, b'(0) ≠ 0 and b(0) = 0. Under these conditions a is invertible and b has a compositional inverse. But it also means that a Riordan square cannot be a Riordan array because a(0) ≠ 0 would contradict b(0) = 0. As we have seen above for the Riordan square the requirement that the constant term and the x term of the series does not vanish is natural and ensures that the Riordan square is invertible as a matrix.

Computationally this means that we could have implemented the Riordan square as follows:

    def RiordanProduct(a, b, dim, expo=false):
        A = toList(a, dim)
        B = A if b == None else toList(b, dim)
        M = matrix(QQ, dim, dim)
        for k in (0..dim-1): M[k, 0] = A[k]
        for k in (1..dim-1):
            for m in (k..dim-1):
                M[m, k] = sum(M[j, k-1]*B[m-j] for j in (k-1..m-1))
        if expo: return ExponentialWeights(M)
        return M 

    def RiordanSquare(a, n, expo=false):
        return RiordanProduct(a, None, n, expo) 

In this setup RiordanSquare(a) = RiordanProduct(a, a), similar to a2 = a*a. The function takes care that the series of a is not expanded twice.

Stieltjes continued fractions

Riordan arrays, continued fractions and Deléham's Δ-operation, which we will discuss below, are closely connected, a fact that is not obvious. Therefore we briefly discuss continued fractions in the form studied by Jacobi and Stieltjes. For details we refer the reader to section 3.10 in the Digital Library of Mathematical Functions, DLMS.

To use continued fractions (CF) as generating functions (GF) for Riordan squares we first define two functions, StieltjesCF and StieltjesGF. The parameter a represents a sequence and the parameter p is a nonnegative integer with default value 2 (the classical case studied by Jacobi).

def StieltjesCF(a, p=2):
    m = 1
    for k in range(len(a) - 1, -1, -1):
        m = 1 - a[k]*x^p/m
    return 1/m

def StieltjesGF(a, p=2):
    cf = StieltjesCF(a, p)
    return cf.series(x, len(a) - 1).list()

For example the sequence S = [n+k for n in (0...)] gives for k = 3 the continued fraction StieltjesCF(S, p=1) = 1/(1-3x/(1-4x/(1-5x/(1-6x/(1-7x/(...)))))). For k ≥ 0 the coefficients of the corresponding generating series start:

    def s(n): return [n+k for k in (0..10)]
    for n in (0..7): print(StieltjesGF(s(n), p=1))

    1, 0,  0,   0,    0,      0,       0,        0,         0, ...
    1, 1,  3,  15,  105,    945,   10395,   135135,   2027025, ...
    1, 2, 10,  74,  706,   8162,  110410,  1708394,  29389186, ...
    1, 3, 21, 207, 2529,  36243,  591381, 10605087, 201756609, ...
    1, 4, 36, 444, 6636, 114084, 2134116, 41682204, 830608236, ...

The fact that the second row of the output matrix are the double factorials of odd numbers (A001147) was proved by Gauß. A321964 records the array. The squares give another noteworthy example: StieltjesGF([n^2 for n in (1..9)], p=1) returns the Euler numbers A000364.

Jacobi continued fractions

Jacobi continued fractions JacobiCF depend on two given sequences a and b. Jacobi generating functions JacobiGF give the series expansion of the Jacobi continued fractions.

    def JacobiCF(a, b, p=2):
        m = 1
        for k in range(len(a)-1, -1, -1):
            m = 1 - b[k]*x - a[k]*x^p/m
        return 1/m

    def JacobiGF(a, b, p=2):
        cf = JacobiCF(a, b, p)
        l = min(len(a), len(b))
        return cf.series(x, l).list()

    def JacobiSquare(a, p=2):
        cf = JacobiCF(a, a, p)
        return cf.series(x, len(a)).list()

The Stieltjes continued fractions are the special case b = 0 of the Jacobi continued fractions. We will call the case a = b the Jacobi square of a and write JacobiSquare(a) = JacobiGF(a, a).

Some examples of Jacobi squares are recorded in A321960. In particular the Bell numbers appear in row 1 in the array below.

    def s(n): return [n+k for k in (0..10)]
    for n in (0..9): print(JacobiSquare(s(n)))

    1, 0,  0,  0,   0,    0,     0,     0,      0,       0, ...
    1, 1,  2,  5,  15,   52,   203,   877,   4140,   21147, ...
    1, 2,  6, 22,  92,  426,  2146, 11624,  67146,  411142, ...
    1, 3, 12, 57, 303, 1752, 10845, 71139, 491064, 3549333, ...

Other examples are the pattern-avoiding permutations A049774 and series-parallel networks with n labeled edges A006351, which correspond to the Bell numbers in the case p = 1.

    r = [n^2 for n in (1..9)] 
    s = [  n for n in (1..9)]

    JacobiGF(r, s)
    1, 1, 2, 5, 17, 70, 349, 2017, 13358, 99377, 822041, ...

    JacobiSquare(s, p=1)
    1, 2, 8, 52, 472, 5504, 78416, 1320064, 25637824, ...

Deléham's Δ-operation

In this section we show an implementation of Philippe Deléham’s Δ-operation, which maps, similar to the Riordan product, two integer sequences on a lower triangular matrix.

    def DelehamDelta(R, S, dim=0):
        x, y = var('x, y')
        r = toList(R, dim) 
        s = toList(S, dim) 
        dim = min(len(r), len(s))
        g = SR(0)
        for n in range(dim-1, -1, -1):
            g = (r[n]*y + s[n]*x*y)/(SR(1) - g)
        ser = (1/(1 - g)).series(y, dim)
        return [expand(p).list() for p in ser.list()]

Analyzing this function demystifies the Δ-operation. The key observation is that in the body of the for-loop 'g' appears on both sides of the assignment. This effectively computes a continued fraction depending on the two input sequences!

Like with the RiordanSquare function we use the polymorphic interface also for the Δ-operation. The next example shows the three possible ways to call the DelehamDelta function:

    def R(n): return 1 if 2.divides(n) else (n+1)//2
    def S(n): return 0^n

    DelehamDelta(R, S, 8)
    DelehamDelta((1 + x - x^2)/(1 - 2*x^2 + x^4), 1, 8)
    DelehamDelta([1, 1, 1, 2, 1, 3, 1, 4], [1, 0, 0, 0, 0, 0, 0, 0])
    
    1,
    1,   1,
    2,   3,    1,
    5,   9,    5,    1,
    15,  29,   20,   7,   1,
    52,  102,  77,   35,  9,   1,
    203, 392,  302,  157, 54,  11, 1,
    877, 1641, 1235, 683, 277, 77, 13, 1. 

The Deléham transform

One can often see an expression like R = S Δ [1, 0, 0, 0, ...] in the OEIS. The expression looks suspiciously similar to a = b * 1 and suggests to eliminate the right factor from the product. This can be easily done.

Let's define the Deléham transform, which depends only on one sequence S and returns a sequence, not a triangle like the Deléham product.

    def DelehamTransform(S, dim):
        x, y = var('x, y')
        s = toList(S, dim)
        g = SR(0)
        for n in range(dim-1, 0, -1):
            g = s[n]*y/(SR(1) - g)
        g = ((s[0] + x)*y)/(SR(1) - g)
        ser = (1/(1 - g)).series(y, dim)
        return [expand(p).list()[0] for p in ser.list()]

We will call S the Deléham generating function of the DelehamTransform(S). For example the Deléham generating function of the Bell numbers is given by the sequence a(n) = 1 for even n and (n+1)/2 for odd n, as we expect from the example in last section. The point is that DelehamTransform(S) returns the first column of the triangle DelehamDelta(S, 1). But this column is the backbone of the Riordan square. This means that there is the identity

    S Δ [1, 0, 0, 0, ...] = RiordanSquare(DelehamTransform(S)).

There are more than 600 sequences in the OEIS defined by Deléham's Δ-operation, many of them in the form on the left hand side of this equation. This identity shows that all these cases are Riordan squares!

The Riordan square of the Bell numbers

Paul Barry, who probably added more Riordan arrays to the OEIS than anyone else, contributed the next example, the Riordan square of the Bell numbers. It starts from 1, 2, 3, ... by first applying the JacobiSquare and then the RiordanSquare transform. The result is A154380.

    S = [n for n in (1..9)]
    RiordanSquare(JacobiSquare(S))

    1;
    1,   1;
    2,   3,    1;
    5,   9,    5,    1;
    15,  29,   20,   7,   1;
    52,  102,  77,   35,  9,   1;
    203, 392,  302,  157, 54,  11, 1;
    877, 1641, 1235, 683, 277, 77, 13, 1;

Compare this triangle with the last one given above: it is the same. Thus in this example the Δ-transform coincides with the JacobiSquare continued fraction transform resulting in the given factorization of the Δ-operation.

Paul Barry comments in this entry (which deserves the rating 'nice') that the Deléham triangle [r0, r1, r2, ...] Δ [s0, s1, s2, ...] has as generating function the bivariate continued fraction

    cf(x, y) = 1/(1 - (r_0*x + s_0*x*y)
                    /(1 - (r_1*x + s_1*x*y)
                        /(1 - (r_2*x + s_2*x*y)
                            /(1 - ...

Observe that this is exactly how we implemented the function DelehamDelta above.

Rational generators for the exponential case

Up to now we assumed our transformations to start from integer sequences. This is an arbitrary limitation since we have in fact only used operations common to all rings. But now we will work with exponential Riordan squares.

For example what happens if we start the exponential RS transformation with the rational sequence 1, 1, 1/2, 1/6, 1/24, 1/120, 1/720, 1/5040, ...? We still get a nice integer triangle! If we remember the generating function of this sequence we can call:

    RiordanSquare(exp(x), 8, true)

    1;
    1,    1;
    1,    3,    1;
    1,    7,    6,    1;
    1,   15,   25,   10,    1;
    1,   31,   90,   65,   15,    1;
    1,   63,  301,  350,  140,   21,    1;
    1,  127,  966, 1701, 1050,  266,   28,    1;

Wow, this is the triangle of Stirling numbers of the second kind, A008277. However, we must add, with the wrong offset. As we have said all our Riordan square triangles are (0,0) -based whereas A008277 is (1,1)-based. But let's not be distracted by this right now and look at another example.

If we start with the rational sequence 0, 1, 1/2, 1/3, 1/4, 1/5, 1/6, 1/7, ... (note that it is not verboten to start the generating sequence with 0, we only cannot expect an invertible triangle in return), then we get:

    RiordanSquare(-ln(1 - x), 8, expo=true)

      0;
      1,    0;
      1,    2,    0;
      2,    6,    3,    0;
      6,   22,   18,    4,    0;
     24,  100,  105,   40,    5,    0;
    120,  548,  675,  340,   75,    6,  0;
    720, 3528, 4872, 2940,  875,  126,  7,  0;

Again this is a well known triangle, A028421. This example has a counterpart only recently studied by W. Lang in A321331.

    RiordanSquare(exp(x) - 1, 8, expo=true)

    0;
    1,    0;
    1,    2,    0;
    1,    6,    3,    0;
    1,   14,   18,    4,   0;
    1,   30,   75,   40,   5,   0;
    1,   62,  270,  260,  75,   6, 0;
    1,  126,  903, 1400, 700, 126, 7, 0;

The main diagonal with the 0's is equivalent to a shift in the index based representation. The last example is defined as A321331(n,k) = (k+1)*Stirling2(n+1, k+1) and the one above A028421(n, k) = (k+1)*Stirling1(n+1, k+1).

Another exponential Riordan square

Surprising few exponential Riordan squares are in the OEIS, less than ten. Among these there is a little known sequence given by Werner Schulte in A265649, a nice family of polynomials. The generating rational sequence is binomial(2n,n)/(2^n) = (2n-1)!!/n!. If you prefer to work with generating series you can call alternatively RiordanSquare((1 - 2*x)^(-1/2), 8, expo=true).

   RiordanSquare([1, 1, 3/2, 5/2, 35/8, 63/8, 231/16, 429/16], expo=true)

        1;
        1,      1;
        3,      5,      1;
       15,     33,     12,      1;
      105,    279,    141,     22,     1;
      945,   2895,   1830,    405,    35,    1;
    10395,  35685,  26685,   7500,   930,   51,  1;
   135135, 509985, 435960, 146685, 23310, 1848, 70,  1;

The inverse of this matrix gives rise to a family of orthogonal polynomials which are now documented in A321966. Phenomenological they resemble very much the Eulerian polynomials. We conjecture that the polynomials have only negative and simple real roots. They deserve a closer study.

OrthogonalPolynomials2
    p(0, x) = 1;
    p(1, x) = x + 1;
    p(2, x) = x^2 + 5*x + 2;
    p(3, x) = x^3 + 12*x^2 + 27*x + 6;
    p(4, x) = x^4 + 22*x^3 + 123*x^2 + 168*x + 24;
    p(5, x) = x^5 + 35*x^4 + 365*x^3 + 1275*x^2 + 1200*x + 120;
    p(6, x) = x^6 + 51*x^5 + 855*x^4 + 5655*x^3 + 13950*x^2 + 9720*x+720;

We might generalize the above considerations by looking at a sequence of polynomial sequences (or a sequence of triangles if you prefer). The recipe is simple, just set m = 1, 2, 3,... in

    RiordanSquare((1 - m*x)^(-1/m), expo=true).inverse().

Case m = 1 gives the Laguerre polynomials in their unsigned form, m = 2 is the case we just discussed and the case m = 3 is:

    p(0, x) = 1;
    p(1, x) = x + 1;
    p(2, x) = x^2 + 6*x + 2;
    p(3, x) = x^3 + 15*x^2 + 38*x + 6;
    p(4, x) = x^4 + 28*x^3 + 188*x^2 + 272*x + 24;
    p(5, x) = x^5 + 45*x^4 + 580*x^3 + 2340*x^2 + 2200*x + 120;
    p(6, x) = x^6 + 66*x^5 + 1390*x^4 + 11040*x^3 + 30280*x^2+19920*x+720;

Production matrices

How can we prove that the polynomials are orthogonal? The classical idea goes back to Favard and Stieltjes who showed (in different ways) that a sequence of polynomials satisfying a suitable 3-term recurrence relation is a sequence of orthogonal polynomials. Here we consider the general case of m -orthogonal polynomials, polynomials which can be computed by a (m + 1)-recurrence. For this purpose we introduce the production matrix of a Riordan square R which exposes the coefficients of the terms in the recurrence.

We need some notation. If M is an invertible matrix M with inverse M(-1) which has top row t then [M−t] is M with row t removed. If W is a matrix and s a row then [W+s] is W with s added as top row. With this notation we define P(R) = [(R(-1) * [R − t]) + e] where t is the top row of R and e is the row [1, 0, 0,...].

Note that usually the production matrix of R is defined as a lower Hessenberg matrix. We define it as a regular lower triangular matrix. Essentially this just adds a 1 at the top of the triangle. This has the advantage that production matrices can be entered in the OEIS as first class triangles and do not have the form of so-called 'funny tables' with cut-off heads, often hidden in the example sections of other sequences.

    def ProductionMatrix(g, dim, expo=false):
        R = RiordanSquare(g, dim+1, expo)
        I = R.inverse()
        r = matrix([R.row(n)[0:dim] for n in (1..dim)])
        i = matrix([I.row(n)[0:dim] for n in (0..dim-1)])
        P = i * r
        R = matrix.identity(dim)
        for n in (0..dim-2): 
            for k in (0..n):
                R[n+1, k] = P[n, k]
        return R

Continuing our example from above we show the production matrices of the exponential Riordan square (1 - m*x)^(-1/m) for the cases m = 3 and m = 4.

   ProductionMatrix((1-3*x)^(-1/3), 8, expo=true)

    [  1   0   0   0   0   0   0   0]
    [  1   1   0   0   0   0   0   0]
    [  3   5   1   0   0   0   0   0]
    [  6  18   9   1   0   0   0   0]
    [  6  42  45  13   1   0   0   0]
    [  0  48 132  84  17   1   0   0]
    [  0   0 180 300 135  21   1   0]
    [  0   0   0 480 570 198  25   1]

    ProductionMatrix((1-4*x)^(-1/4), 8, expo=true)

    [   1    0    0    0    0    0    0    0]
    [   1    1    0    0    0    0    0    0]
    [   4    6    1    0    0    0    0    0]
    [  12   28   11    1    0    0    0    0]
    [  24   96   72   16    1    0    0    0]
    [  24  216  312  136   21    1    0    0]
    [   0  240  840  720  220   26    1    0]
    [   0    0 1080 2280 1380  324   31    1]

The case m = 3 is recorded in A322944. We refer the reader to this entry to see that the rows of the production matrix give the coefficients of the recurrence of the polynomials. A cool result.

Riordan squares in the OEIS

Some Riordan squares which can be found in the OEIS are given in the table below. See also A321620.

Riordan Squares
Sequence a RS(a) gf(a)
Catalan A039599 (1 - sqrt(1 - 4x))/(2x)
Motzkin A321621 (1 - x - sqrt(1 - 2x - 3x^2))/(2*x^2)
Fine A321622 1 + (1 - sqrt(1 - 4x))/(3 - sqrt(1 - 4x))
large Schröder A321623 (1 - x - sqrt(1 - 6x + x^2))/(2x)
little Schröder A172094 (1 + x - sqrt(1 - 6x + x^2))/(4x)
Lucas A321624 1 + x(1 + 2x)/(1 - x - x^2)
orbital numbers A321625 (1 + x/(1 - 4x^2))/sqrt(1 - 4x^2)
ternary trees |A109956| u = sqrt(x3/4); sin(arcsin(3u)/3)/u
central trinomial A116392 1/sqrt(1 - 2x - 3x^2)
1, 2, 3, 4, ... A210753 1/(1-x)^2
powers of 2 A038208 1/(1 - 2*x)
the all 1's seq. A007318 1/(1 - x)
Fibonacci A063967 1/(1 - x - x^2)
Tribonacci A187889 1/(1 - x - x^2 - x^3)
(2*n-1)!! A321627 StieltjesGF(n ↦ n+1, p=1)
Euler A321630 StieltjesGF(n ↦ (n+1)^2, p=1)
involutions A321629 JacobiGF(n ↦ n+1, n ↦ 1, p=2)
Bell A154380 JacobiSquare(n ↦ n+1, p=2)
Exponential RS
Laguerre |A021009| 1/(1 - x)
Stirling2 A008277 exp(x)
Narumi A028421 -ln(1 - x)
W. Schulte A265649 (1 - 2*x)^(-1/2)
Orthogonal polynomials
new A321966 inverse matrix (1-2*x)^(-1/2)
new A322944 inverse matrix (1-3*x)^(-1/3)
Jacobsthal A322942 (2*x^2-1)/((x + 1)*(2*x - 1))

We also added A322942, the coefficients of orthogonal polynomials J(n, x) which generalize the Jacobsthal numbers with J(0) = 1. Also noteworthy is the sequence 2n J(n, 1/2) reported in A323232.

There are more Riordan square arrays in the OEIS. For example what W. Lang calls Pascal (1,n) triangles (A095660, A095666, A096940) are Riordan squares.

Searching for "Riordan array" in the OEIS gives currently over 1300 hits and there is a list of Riordan arrays in the OEIS Wiki compiled by Ralf Stephan which seeks to give a systematic overview.

Hints to the literature

As a starter we recommend Renzo Sprugnoli's classic "Riordan arrays and combinatorial sums" from the 1990s. We further recommend Sprugnoli's "An Introduction to Mathematical Methods in Combinatorics" which you can find in the Internet. It has a chapter on Riordan arrays.

We allude to the recent (2017) book "Riordan Arrays: A Primer" from Paul Barry.

There are two YouTube videos from "Experimental Mathematics" if you like this presentation form. Search for "Riordan Arrays and Their Applications in Combinatorics", part 1 and 2.

Appendix: An interface for the Riordan and Deléham functions

We compiled the functions given above in a SageMath/Jupyter notebook which you can download from GitHub. Error reports and pull requests are welcome.

    def TaylorList(f, n):
        t = f.taylor(x, 0, n+1).list()
        return (t + [0]*(n - len(t)))[0:n]

    def toList(R, dim):  
        if isinstance(R, Integer): R = SR(R)
        if isinstance(R, type(SR())):
            if dim == 0: raise ValueError('dim must not be zero.') 
            return TaylorList(SR(R), dim)
        if callable(R):
            if dim == 0: raise ValueError('dim must not be zero.') 
            return [R(n) for n in (0..dim-1)]
        return R