Multivariate Pade approximation, recognizing rational function

edit

(Question edited)

Suppose that I have a taylor polynomial approximation of the taylor series of a rational function $f(x,y) =p(x,y)/q(x,y)$

that I want to recover.

1) If $f(x,y) \in Q[y][[x]]$ then the following approach seems to work: convert to $Q[y][[x]]$ which magically gets embedded into $Q(y)[[x]]$, and then standard pade approximation works.

S.<y> = QQ[]
R.<x> = S[]
U.<x> = PowerSeriesRing(S,default_prec=7)
p = x+y
q = (1-x)*(1-x*y)
fcorrect=U(p)/U(q)
fperturbed = fcorrect + U(x^7*y^7)
print(fperturbed)
print()
print(fperturbed.pade(3,3))

output:

y + (y^2 + y + 1)*x + (y^3 + y^2 + 2*y + 1)*x^2 + (y^4 + y^3 + 2*y^2 + 2*y + 1)*x^3 + (y^5 + y^4 + 2*y^3 + 2*y^2 + 
2*y + 1)*x^4 + (y^6 + y^5 + 2*y^4 + 2*y^3 + 2*y^2 + 2*y + 1)*x^5 + (y^7 + y^6 + 2*y^5 + 2*y^4 + 2*y^3 + 2*y^2 + 
2*y + 1)*x^6 + O(x^7)

(x + y)/(y*x^2 + (-y - 1)*x + 1)

2) On the other hand, for rational power series like $1/((1-x)(1-y))$ which, as power series in x, have coefficients that are rational functions in y rather than polynomials in y, I had to manually convert to $Q(y)[[x]]$:

b=(1-x)*(1-y)
V.<x> = PowerSeriesRing(S.fraction_field(),default_prec=9)
T.<x,y> = PowerSeriesRing(QQ,'x,y',default_prec=9)
W.<y> = PowerSeriesRing(QQ,'y',default_prec=9)

hcorrect = 1/T(b)
hperturbed = hcorrect + T(x^7*y^7)
print(hperturbed)
print()
hp = hperturbed.polynomial()
W.<y> = PowerSeriesRing(QQ,'y')

hconverted = add(V(x^k*W(hp.coefficient({x:k})).pade(2,2)) for k in range(9))
print(hconverted)
print()
print(hconverted.pade(2,2))

gives

1 + x + y + x^2 + x*y + y^2 + x^3 + x^2*y + x*y^2 + y^3 + x^4 + x^3*y + x^2*y^2 + x*y^3 + y^4 + x^5 + x^4*y + 
x^3*y^2 + x^2*y^3 + x*y^4 + y^5 + x^6 + x^5*y + x^4*y^2 + x^3*y^3 + x^2*y^4 + x*y^5 + y^6 + x^7 + x^6*y + 
x^5*y^2 + x^4*y^3 + x^3*y^4 + x^2*y^5 + x*y^6 + y^7 + x^8 + x^7*y + x^6*y^2 + x^5*y^3 + x^4*y^4 + x^3*y^5 + 
x^2*y^6 + x*y^7 + y^8 + O(x, y)^9

-1/(y - 1) + (-1/(y - 1))*x + (-1/(y - 1))*x^2 + (-1/(y - 1))*x^3 + (-1/(y - 1))*x^4 + (-1/(y - 1))*x^5 + (y^2 + y + 1)*x^6 + 
(y + 1)*x^7 + x^8

(1/(y - 1))/(x - 1)

The first case 1) is easy enough, and I appreciate the automagic conversion to $Q(y)[[x]]$. It is somewhat slow when the method is extended to several variables, though.

The second case 2) was very hard to get to work; is there an easier way? How do I speed up the calculations?