How to improve ( mathematica ) performance when dealing with recursive functions?

Background. I want to print a table of convergents for 31^(1/2). I made the following recursive definition of the table. ( Exchange 31^(1/2) with the golden ratio and the table below would contain the Fibonacci series ).

 cf := ContinuedFraction
 tf := TableForm
 p[-1] = 0; p[0] = 1; q[-1] = 1; q[0] = 0;
 a[k_] := cf[Sqrt[31], k][[k]]
 p[k_] := a[k]*p[k - 1] + p[k - 2]
 q[k_] := a[k]*q[k - 1] + q[k - 2]
 s[n_] := Timing[Table[{k, a[k], p[k], q[k]}, {k, 8, 8 n, 8}]] // tf

Timing increases exponentially fast. I had to alt+. (abort) for s[4].

Question: How to improve ( mathematica ) performance when dealing with recursive function开发者_运维技巧s?

From a quick (not thorough, to admit it) look at your code, it looks like both p and q are defined recursively in terms of two previous values. This means that to calculate the nth value of p, ~2^n evaluations are needed (every step doubles the number). So yes, complexity is exponential, regardless of Mathematica or any other language being used.

If you insist on using a recursive formulation of the problem (e.g. for simplicity), then the simplest way to reduce the performance penalty is using memoization, i.e. doing something like

p[k_] := p[k] = a[k]*p[k - 1] + p[k - 2]

Don't forget to Clear[p] before any redefinition.

In short, memoization means having the function remember the computation result for each input, so subsequent evaluations are faster. It's very likely faster, but more complicated, to compute two values (p_(n+1) and p_(n)) from two previous values (p_(n) and p_(n-1)), then the complexity will be linear instead of exponential.

_{I hope this helps. I don't have Mathematica here to test right now.}

Here is a small further refinement. Since this is a quadratic irrational you can also compute the a[k] coefficients more directly.

In[499]:= Clear[a, p, q, cf]
cf = ContinuedFraction[Sqrt[31]];
cf2len = Length[cf[[2]]];
a[1] = cf[[1]];
a[k_] := cf[[2, Mod[k - 1, cf2len, 1]]]
p[-1] = 0; p[0] = 1; q[-1] = 1; q[0] = 0;
p[k_] := p[k] = a[k]*p[k - 1] + p[k - 2]
q[k_] := q[k] = a[k]*q[k - 1] + q[k - 2]
s[n_] := Timing[Table[{k, a[k], p[k], q[k]}, {k, 8, 8 n, 8}];]

In[508]:= s[1000]

Out[508]= {0.12, Null}

In[509]:= Clear[a, p, q, cf]
cf := ContinuedFraction
p[-1] = 0; p[0] = 1; q[-1] = 1; q[0] = 0;
a[k_] := a[k] = cf[Sqrt[31], k][[k]]
p[k_] := p[k] = a[k]*p[k - 1] + p[k - 2]
q[k_] := q[k] = a[k]*q[k - 1] + q[k - 2]
s[n_] := Timing[Table[{k, a[k], p[k], q[k]}, {k, 8, 8 n, 8}];]

In[516]:= s[1000]

Out[516]= {6.08, Null}

Also you can get a[k] in closed form, though it is not terribly pretty.

In[586]:= Clear[a];
asoln[k_] = 
 FullSimplify[
  a[k] /. First[
    RSolve[Join[
      Table[a[k] == cf[[2, Mod[k - 1, cf2len, 1]]], {k, 
        cf2len}], {a[k] == a[k - 8]}], a[k], k]], Assumptions -> k > 0]

Out[587]= (1/(8*Sqrt[2]))*(4*(Cos[(k*Pi)/4] + Sin[(k*Pi)/4])*
        (-2*Sqrt[2] + (5 + 2*Sqrt[2])*Sin[(k*Pi)/2]) + 
      Sqrt[2]*(25 - 9*Cos[k*Pi] + 26*Sin[(k*Pi)/2] - 9*I*Sin[k*Pi]))

Offhand I do not know whether this can be used to get a direct solution for p[k] and q[k]. RSolve seems unable to do that.

--- edit ---

As others have mentioned, it can be cleaner to just build the list from first to last. Here is the handling of p[k], using memoization as above vs NestList.

Clear[a, p, q, cf]
cf = ContinuedFraction[Sqrt[31]];
cf2len = Length[cf[[2]]];
a[1] = cf[[1]];
a[k_] := cf[[2, Mod[k - 1, cf2len, 1]]]

p[-1] = 0; p[0] = 1;
p[k_] := p[k] = a[k]*p[k - 1] + p[k - 2]
s[n_] := Timing[Table[p[k], {k, n}];]

In[10]:= s[100000]

Out[10]= {1.64, Null}

In[153]:= s2[n_] := Timing[ll = Module[{k = 0},
      NestList[(k++; {#[[2]], a[k]*#[[2]] + #[[1]]}) &, {0, 1}, 
       n]][[All, 2]];]

In[154]:= s2[100000]

Out[154]= {0.78, Null}

In addition to being somewhat faster, this second approach does not keep a large number of definitions around. And you do not really need them in order to generate more elements, because this iteration can be resumed using a pair from the last elements (make sure they start at 0 and 1 modulo 8).

I will mention that one can obtain a closed form for p[k]. I found it convenient to break the solution into 8 (that is, cf2len) pieces and link them via recurrences. The reasoning behind the scenes comes from basic generating function manipulation. I did some slightly special handling of one equation and one initial condition to finesse the fact that a[1] is not part of the repeating sequence.

In[194]:= func = Array[f, cf2len];
args = Through[func[n]];
firsteqns = {f[2][n] == a[2]*f[1][n] + f[cf2len][n - 1], 
   f[1][n] == a[9]*f[cf2len][n - 1] + f[cf2len - 1][n - 1]};
resteqns = 
  Table[f[j][n] == a[j]*f[j - 1][n] + f[j - 2][n], {j, 3, cf2len}];
inits = {f[8][0] == 1, f[1][1] == 5};
eqns = Join[firsteqns, resteqns, inits];

In[200]:= 
soln = FullSimplify[args /. First[RSolve[eqns, args, n]], 
   Assumptions -> n > 0];

In[201]:= FullSimplify[Table[soln, {n, 1, 3}]]

Out[201]= {{5, 6, 11, 39, 206, 657, 863, 1520}, {16063, 17583, 33646, 
  118521, 626251, 1997274, 2623525, 4620799}, {48831515, 53452314, 
  102283829, 360303801, 1903802834, 6071712303, 7975515137, 
  14047227440}}

Quick check:

In[167]:= s2[16]; ll

Out[167]= {1, 5, 6, 11, 39, 206, 657, 863, 1520, 16063, 17583, 33646, \
118521, 626251, 1997274, 2623525, 4620799}

We can now define a function from this.

In[165]:= 
p2[k_Integer] := soln[[Mod[k, cf2len, 1]]] /. n -> Ceiling[k/cf2len]

In[166]:= Simplify[p2[4]]

Out[166]= 39

I do not claim that this is particularly useful, just wanted to see if I could actually get something to work.

--- end edit ---

Daniel Lichtblau