Using Fold to calculate the result of linear recur

Question

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 n^th 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 2¹⁰ - 2¹³ 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, {z_i}. Part of creating the approximant is to generate a set of coefficients, a_i, 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.

Answer 1

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}

Answer 2

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[]).

Answer 3

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)

Answer 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