How to compute Riordan arrays?

savecancel

def riordan_array(d, h, n, exp=false):
   """
   The function computes the Riordan array of the
   formal power series d and h as a lower triangular
   matrix of dimension n.

   If the parameter 'exp' is true the corresponding
   exponential Riordan array is computed.
   """

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

   td = taylor_list(d, n)
   th = taylor_list(h, n)
   M = matrix(QQ, n, n)

   for k in (0..n-1): M[k, 0] = td[k]

   for k in (1..n-1):
       for m in (k..n-1):
           M[m, k] = add(M[j, k-1]*th[m-j] for j in (k-1..m-1))

   if exp:
       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

We give three examples:

riordan_array(1/(1-x), x/(1-x), 7)
[ 1  0  0  0  0  0  0]
[ 1  1  0  0  0  0  0]
[ 1  2  1  0  0  0  0]
[ 1  3  3  1  0  0  0]
[ 1  4  6  4  1  0  0]
[ 1  5 10 10  5  1  0]
[ 1  6 15 20 15  6  1]

riordan_array(1, log((1-x)/(1-2*x)), 6, true)
[ 1   0   0   0   0  0]
[ 0   1   0   0   0  0]
[ 0   3   1   0   0  0]
[ 0  14   9   1   0  0]
[ 0  90  83  18   1  0]
[ 0 744 870 275  30  1]

# Fibonacci-Catalan triangle
riordan_array(1/(1-x-x^2), (1-sqrt(1-4*x))/2, 7)
[ 1  0  0  0  0  0  0]
[ 1  1  0  0  0  0  0]
[ 2  2  1  0  0  0  0]
[ 3  5  3  1  0  0  0]
[ 5 12  9  4  1  0  0]
[ 8 31 26 14  5  1  0]
[13 85 77 46 20  6  1]