Fastest way to simplify rational functions

Question

I am using Series to approximate function of two variables:

Series[f[x,y],{x,0,m}]

the function is a complicated sum of hypergeometric functions, but my question is more generic than that. The upshot is that I end up with a sum of rational functions in y. Schematically

$$f(x,y)=\sum_{i=0}^{m} \frac{p_i(y)}{q_i(y)} x^i+\mathcal{O(x^{m+1})}\sim \frac{p(x,y)}{Q(y)}$$

where $p_i,q_i$ are polynomials of $y$ only. I know what the all the denominator divide some polynomial $Q(y)$m i.e. the roots of $q_i(y)$ are in the roots, $a_j$, of $Q(y)$ (with multiplicity). I need to know the result as a rational function $p(x,y)/Q(y)$ in order to use some custom functions that only works on polynomials, i.e. they use e.g. CoefficientList and other polynomial-related functions.

For illustration let us take the function as a simple hypergeometric function:

(* f[x,y] is more complicated than that. For illustration only.*)
m = 10;
f[x_,y_]   := Hypergeometric2F1[y, y + 1, y - 3, x/(x^2 + 1)];
lista       = {1,2,3};
denominator = Product[y-a,{a,lista}];
numerator   = Series[denominator * f[x,y],{x,0,m}]

In principle, numerator should be nice polynomial in both $x,y$. In practice however, what I get are sums complicated rational functions that Mathematica does not simplify. For low m this is not an issue, and I simply do

Collect[numerator//Normal,x,Simplify]

And I get the result I want. If I expand the series to high orders however, this quickly becomes very long and slows down the whole script. Is there a best practice for simplifying rational functions like this? I tried the following:

Collect[numerator//Normal,x,Simplify]
Collect[numerator//Normal,x,Together]
Simplify[numerator[[3]]].Table[x^i, {i, 0, m}]
Together[numerator[[3]]].Table[x^i, {i, 0, m}]

They all take the same order of magnitude in time to process at large m.

On my laptop, for this example function and m = 100, the series is found in about 10 seconds and all the above each take between 50 and 60 seconds to get the desired result.

What is the best practice and faster way in this case?

Answer 1

It can be done somewhat faster by interpolating the polynomials in y. I'll illustrate with an expansion to degree 64.

m = 64;
f[x_, y_] := Hypergeometric2F1[y, y + 1, y - 3, x/(x^2 + 1)];
lista = {1, 2, 3};
denominator = Product[y - a, {a, lista}];
numerator = Series[denominator*f[x, y], {x, 0, m}];

(1) Just use Together:

Timing[t1 = Together[Normal@numerator];]
LeafCount[t1]

(* Out[455]= {7.48, Null}

Out[456]= 18352 *)

(2) Use Together on individual "coefficients' in y. Also make sure results are equivalent.

Timing[
 t2 = Together[numerator[[3]]].Table[x^i, {i, 0, m}];]
LeafCount[t2]
Together[t2 - t1]

(* Out[475]= {4.74, Null}

Out[476]= 11472

Out[477]= 0 *)

(3) Interpolate in y for each coefficient to clear denominators. For this I decided to use integers starting 10 larger than m and to use around 3/2 m of them. Not sure this many are needed though, or if this is in general an adequate range to use.

Timing[
 t3 = numerator[[3]];
 t3b = Transpose[Table[t3, {y, 10 + m, 10 + Ceiling[5/2*m]}]];
 t3c = Map[
   Together[
     InterpolatingPolynomial[
      Transpose[{Range[10 + m, 10 + Ceiling[5/2*m]], #}], y]] &, 
   t3b];
 t3d = t3c.Table[x^i, {i, 0, m}];]
LeafCount[t3d]
Together[t3d - t2]

(* Out[480]= {1.692, Null}

Out[481]= 11679

Out[482]= 0 *)

So it's a modest improvement, factor of 3 or so at this range for m.