# Speeding up ReplaceRepeated while truncating to desired order

I need to program in an algorithm that recursively makes algebraic replacements which leads to an utterly complicated algebraic function of $x$, but whose final result is only needed at fixed order in the Laurent series in $x$. I wish to use this fact to make the algorithm go faster. Here is sample problem that mimics the realistic problem:

\begin{align}f_n(x) &= \sum_{j=1}^4\frac{1}{6+j x}n\left(\frac{1}{x}+j \,f_{n-1}(x)\right); \qquad f_0(x) = {\tt f[0]}\\ F(x) &= \sum_{n=1}^6\frac{1}{2+n x}\left(\frac{n}{x}+f(n)\right) \end{align}

The algorithm is to repeatedly insert the first equation into the second one until all instances of $f_n(x)$ turn into $f_0={\tt f[0]}$. All that is needed in the end is the Laurent expansion out to $\mathcal{O}(x^0)$ (just need the $1/x$ and constant terms). nota bene: The sums are there only to lengthen the expression to make Mathematica work harder.

Here is a naive implementation (with the result on my machine):

(*Implementation 1*)
f[n_] /; n > 0 -> Sum[1/(6 + j x) n (1/x + j f[n - 1]), {j, 1, 4}];
Sum[1/(2 + n x) (n/x + f[n]), {n, 1, 6}] //. %;
Timing[Series[%, {x, 0, 0}]]


Performing the actual recursion in the second line is lightening fast. But the final step Series function is painfully slow (4.109 seconds), since the $x$ in the denominators forces Mathematica to appeal to its analytical routines to expand things like $\frac{1}{1+x}\approx 1 -x$.

Since all that I need is the first two terms of the Laurent series, I went ahead and, by hand, expanded the fractions $\frac{1}{6+jx}\approx(\frac{1}{6}-\frac{j x}{36})$ and $\frac{1}{2+nx}\approx(\frac{1}{2}-\frac{n x}{4})$, with the intention of partially freeing Mathematica of painful analytics.

(*Implementation 2*)
f[n_] /; n > 0 -> Sum[(1/6 - (j x)/36) n (1/x + j f[n - 1]), {j, 1, 4}];
Sum[(1/2 - (n x)/4) (n/x + f[n]), {n, 1, 6}] //. %;
Timing[Series[%, {x, 0, 0}]]


This is 8 times faster; the more I save Mathematica from doing analytic manipulations, the faster it goes. But I need it to go faster, yet. How can I carry out the alogorithm to minimize the extent Mathematica has to do computationally costly analytic work?

Here's my failed attempt: ask Mathematica drop higher order terms at each stage in the recursion by adding an O[x] to the end of the line:

(*Implementation 3 [FAILURE]*)
f[n_] /; n > 0 -> Sum[(1/6 - (j x)/36) n (1/x + j f[n - 1]), {j, 1, 4}];
Timing[Sum[(1/2 - (n x)/4) (n/x + f[n]), {n, 1, 6}] + O[x] //. %]


Very fast, but the wrong answer. It gets the wrong answer because it drops terms like x*f[n] as it thinks f[n] is order O[x]^0 when in reality there is a 1/x part. What can I do to speed this up?

Truncate early and often.

f[0] = f0;
f[n_] /; n > 0 :=
f[n] = Sum[
Series[1/(6 + j x) n (1/x + j f[n - 1]), {x, 0, 0}], {j, 1, 4}];

Timing[
res = Series[
Sum[Series[1/(2 + n x) (n/x + f[n]), {x, 0, 0}], {n, 1, 6}], {x, 0,
0}]]

(* Out[140]= {0.012000, SeriesData[x, 0, {
Rational[1016635, 162], Rational[5, 486] (-3351289 + 835701 f0)}, -1,
1, 1]} *)

• Note to self: the key is not to use ReplaceRepeated but to assign DownValues – QuantumDot Oct 1 '14 at 9:50