4949.sagews

from sage.symbolic.expression_conversions import polynomial
var('x','y')

#f is the function g is the powerseries expansion wrt variable, of the desired order
f(x,y)=(1+x+y*x^2)/(1-x-x^2-y*x^3)
g(x,y)= f(x,y).series(x,10)

# this gives you the coefficients of the monomials in a vector
# note that anything evaluated to 0 is omitted so this vector is a little useless... 
g(x,y).polynomial(SR).coefficients()

[4, 51, 1, 130, 22, 89, 71, 9, 55, 38, 3, 34, 20, 1, 21, 10, 13, 5, 8, 2, 5, 1, 3, 2, 1]

# you can look at the coefficient for a specific vector (n,d) [where we have x^n *y^d as in the paper]
# here I compute the coefficient of x^3*y which should be 8 [see page 29 of Tia's thesis]
g(x,y).polynomial(SR).coefficient([1,0])

# a loop to print things... 
n=0
d=0
while n<11:
     print (n,d), "coefficient is", g(x,y).polynomial(SR).coefficient([n,d])
     n=n+1

(0, 0) coefficient is 1 (1, 0) coefficient is 2 (2, 0) coefficient is 3 (3, 0) coefficient is 5 (4, 0) coefficient is 8 (5, 0) coefficient is 13 (6, 0) coefficient is 21 (7, 0) coefficient is 34 (8, 0) coefficient is 55 (9, 0) coefficient is 89 (10, 0) coefficient is 0

# you can look at g itself if you want to hand-verify the coefficients 
print g(x,y)

2*(32*(y + 2)*y + 16*(y + 4)*y + 4*((y + 4)*y + 8*y + 16)*y + 4*(2*(y + 2)*y + (y + 4)*y + 8*y + 16)*y + (8*(y + 2)*y + 4*(y + 4)*y + ((y + 4)*y + 8*y + 16)*y + 32*y + 64)*y + (8*(y + 2)*y + 4*(y + 4)*y + ((y + 4)*y + 8*y + 16)*y + (2*(y + 2)*y + (y + 4)*y + 8*y + 16)*y + 32*y + 64)*y + 128*y + 256)*x^9 + (32*(y + 2)*y + 16*(y + 4)*y + 4*((y + 4)*y + 8*y + 16)*y + 4*(2*(y + 2)*y + (y + 4)*y + 8*y + 16)*y + (8*(y + 2)*y + 4*(y + 4)*y + ((y + 4)*y + 8*y + 16)*y + 32*y + 64)*y + 128*y + 256)*x^8 + 2*(8*(y + 2)*y + 4*(y + 4)*y + ((y + 4)*y + 8*y + 16)*y + (2*(y + 2)*y + (y + 4)*y + 8*y + 16)*y + 32*y + 64)*x^7 + (8*(y + 2)*y + 4*(y + 4)*y + ((y + 4)*y + 8*y + 16)*y + 32*y + 64)*x^6 + 2*(2*(y + 2)*y + (y + 4)*y + 8*y + 16)*x^5 + ((y + 4)*y + 8*y + 16)*x^4 + 4*(y + 2)*x^3 + (y + 4)*x^2 + 2*x + 1