3

I have a linear recurrence problem where the next element relies on more than just the prior value, e.g. the Fibonacci sequence. One method calculating the nth element is to define it via a function call, e.g.

Fibonacci[0] = 0; Fibonacci[1] = 1;
Fibonacci[n_Integer?Positive] := Fibonacci[n] + Fibonacci[n - 1]

and for the sequence I'm working with, that is exactly what I do. (The definition is inside of a Module so I don't pollute Global`.) However, I am going to be using this with 210 - 213 points, so I'm concerned about the extra overhead when I just need the last term and none of the prior elements. I'd like to use Fold to do this, but Fold only passes the immediately prior result which means it is not directly useful for a general linear recurrence problem.

I'd like a pair of functions to replace Fold and FoldList that pass a specified number of prior sequence elements to the function, i.e.

In[1] := MultiFoldList[f, {1,2}, {3,4,5}] (* for lack of a better name *)
Out[1]:= {1, 2, f[3,2,1], f[4,f[3,2,1],2], f[5,f[4,f[3,2,1],2],f[3,2,1]]}

I had something that did this, but I closed the notebook prior to saving it. So, if I rewrite it on my own, I'll post it.

Edit: as to why I am not using RSolve or MatrixPower to solve this. My specific problem is I'm performing an n-point Pade approximant to analytically continue a function I only know at a set number of points on the imaginary axis, {zi}. Part of creating the approximant is to generate a set of coefficients, ai, which is another recurrence relation, that are then fed into the final relationship

A[n+1]== A[n] + (z - z[[n]]) a[[n+1]] A[n-1]

which is not amenable to either RSolve nor MatrixPower, at least that I can see.

2
  • 1
    Is the built in Pade approximant no good in the case that you want?
    – Simon
    Mar 14, 2011 at 4:50
  • 1
    The problem is my function is numerical, i.e. I only know it at certain points, and the built-in function approximates known functions. So, I can't give it a list and have it give me a function back. Nor, do the other built-in interpolating functions work as they don't do rational interpolation, and my function will most likely have poles. So, rational interpolation is necessary.
    – rcollyer
    Mar 14, 2011 at 5:01

4 Answers 4

5

Can RecurrenceTable perform this task for you?

Find the 1000th term in a recurrence depending on two previous values:

In[1]:= RecurrenceTable[{a[n] == a[n - 1] + a[n - 2], 
  a[1] == a[2] == 1}, a, 
   {n, {1000}}]

Out[1]= {4346655768693745643568852767504062580256466051737178040248172\
9089536555417949051890403879840079255169295922593080322634775209689623\
2398733224711616429964409065331879382989696499285160037044761377951668\
49228875}

Edit: If your recurrence is defined by a function f[m, n] that doesn't like to get evaluated for non-numeric m and n, then you could use Condition:

In[2]:= f[m_, n_] /; IntegerQ[m] && IntegerQ[n] := m + n

The recurrence table in terms of f:

In[3]:= RecurrenceTable[{a[n] == f[a[n - 1], a[n - 2]], 
  a[1] == a[2] == 1}, a, {n, {1000}}]

Out[3]= {4346655768693745643568852767504062580256466051737178040248172\
9089536555417949051890403879840079255169295922593080322634775209689623\
2398733224711616429964409065331879382989696499285160037044761377951668\
49228875}
3
  • +1, I need to become more familiar with the "recent" additions to the standard functions. It does exactly what I need it to. I had to wrap it in Quiet as the recurrence references external lists and Part just does not like having an undefined variable passed to it ...
    – rcollyer
    Mar 14, 2011 at 4:19
  • Interesting use of condition. I never think to use it quite like that. For the edit alone, I wish I could add another +1.
    – rcollyer
    Mar 14, 2011 at 4:37
  • 2
    (+1) Just a small comment that checking Heads is faster than pattern matching, so f[m_Integer,n_Integer] would be faster than f[m_, n_] /; IntegerQ[m] && IntegerQ[n].
    – Simon
    Mar 14, 2011 at 4:47
4

A multiple foldlist can be useful but it would not be an efficient way to get linear recurrences evaluated for large inputs. A couple of alternatives are to use RSolve or matrix powers times a vector of initial values.

Here are these methods applied to example if nth term equal to n-1 term plus two times n-2 term.

f[n_] =  f[n] /. RSolve[{f[n] == f[n - 1] + 2*f[n - 2], f[1] == 1, f[2] == 1},
  f[n], n][[1]]

Out[67]= 1/3 (-(-1)^n + 2^n)

f2[n_Integer] := Last[MatrixPower[{{0, 1}, {2, 1}}, n - 2].{1, 1}]

{f[11], f2[11]}

Out[79]= {683, 683}

Daniel Lichtblau Wolfram Research

1
  • +1, for suggesting two reasonable approaches that I'd use if my specific case allowed it. Trust me, I really wanted to!
    – rcollyer
    Mar 14, 2011 at 4:39
3

Almost a convoluted joke, but you could use a side-effect of NestWhileList

fibo[n_] := 
  Module[{i = 1, s = 1}, 
   NestWhileList[ s &, 1, (s = Total[{##}]; ++i < n) &, 2]];  

Not too bad performance:

In[153]:= First@Timing@fibo[10000]
Out[153]= 0.235  

By changing the last 2 by any integer you may pass the last k results to your function (in this case Total[]).

0
2

LinearRecurrence and RecurrenceTable are very useful.

For small kernels, the MatrixPower method that Daniel gave is the fastest.

For some problems these may not be applicable, and you may need to roll your own.

I will be using Nest because I believe that is appropriate for this problem, but a similar construct can be used with Fold.

A specific example, the Fibonacci sequence. This may not be the cleanest possible for that, but I believe you will see the utility as I continue.

fib[n_] :=
  First@Nest[{##2, # + #2} & @@ # &, {1, 1}, n - 1]

fib[15]

Fibonacci[15]

Here I use Apply (@@) so that I can address elements with #, #2, etc., rathern than #[[1]] etc. I use SlotSequence to drop the first element from the old list, and Sequence it into the new list at the same time.

If you are going to operate on the entire list at once, then a simple Append[Rest@#, ... may be better. Either method can be easily generalized. For example, a simple linear recurrence implementation is

 lr[a_, b_, n_Integer] := First@Nest[Append[Rest@#, a.#] &, b, n - 1]

 lr[{1,1}, {1,1}, 15]

(the kernel is in reverse order from the built in LinearRecurrence)

1
  • +1, this looks a lot like the method I came up with for Fold directly. And the Append[Rest@#] method is really nice.
    – rcollyer
    Mar 14, 2011 at 12:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.