# How to compute Riordan arrays?

Riordan arrays are an important tool in the study of sequences. See the examples on the OEIS Wiki. With Maple the computation of Riordan arrays is straightforward, here an example of an exponential Riordan array:

Riordan := (d,h,n,k) -> coeftayl(d*h^k,x=0,n)*n!/k!:
A131222 := (n,k) -> Riordan(1,log((1-x)/(1-2*x)),n,k):
seq(print(seq(A131222(n,k),k=0..n)),n=0..5);

1
0,   1
0,   3,   1
0,  14,   9,   1
0,  90,  83,  18,  1
0, 744, 870, 275, 30, 1


It is not clear to me how Riordan arrays are computed efficiently with Sage. Any hints are appreciated.

edit retag close merge delete

The logical next step would be to open an enhancement ticket.

Sort by » oldest newest most voted
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 + *(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]

more

If f is a symbolic function, you can get its Taylor expansion as follows:

sage: f = log((1-x)/(1-2*x))
sage: f.taylor(x,0,5)
31/5*x^5 + 15/4*x^4 + 7/3*x^3 + 3/2*x^2 + x


Then you can get its k^th coefficient as follows:

sage: f.taylor(x,0,5).coefficient(x^3)
7/3


Or:

sage: f.taylor(x,0,5).coefficients()
[[1, 1], [3/2, 2], [7/3, 3], [15/4, 4], [31/5, 5]]
sage: f.taylor(x,0,5).list()
[0, 1, 3/2, 7/3, 15/4, 31/5]

more

Unfortunately, I can not see how your comments address my question. Can you please elaborate, possibly with a runnable function reproducing the above output?