Ask Your Question

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

Comments

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

( 2015-04-17 08:40:26 +0200 )edit

2 Answers

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 + [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]

more

Comments

Very useful. Many thanks ;)

( 2015-04-18 09:40:43 +0200 )edit

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

Comments

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

( 2015-04-17 10:37:55 +0200 )edit

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Stats

Asked: 2015-04-16 22:03:56 +0200

Seen: 839 times

Last updated: Apr 17 '15