How to deal with recursion formula in Mathematica?

Question

In engineering problems, I am always seeing many recursion formula.

For instance, In the book "The NURBS book", I discovered many recursion formula

• Fibonacci $$f(n+2)=f(n+1)+f(n)$$

• deCasteljau algorithm $$\vec{P}_{k,i}(u_0)=(1-u_0)\vec{P}_{k-1,i}(u_0)+u_0\vec{P}_{k-1,i+1}(u_0)$$

• Horner

And so on.

Previously, I know many Mathematica user used the solution as beow to code.

 fib[n_] := fib[n] = fib[n - 1] + fib[n - 2];
 fib[1] = fib[2] = 1;

 fib[200]

 280571172992510140037611932413038677189525

However, when n became larger, like n=300

fib[300]

$RecursionLimit::reclim: Recursion depth of 256 exceeded. >>

The built-in funcion Fibonacci can give the result normally.

Fibonacci[300]

222232244629420445529739893461909967206666939096499764990979600

Question:

Is it possible to have a general(uniform) method(Namely, like Funcion Programing paradigm like Nest,NestList,FoldList) to deal with these recursion formula in Mathematica?

Summarize

Firstly, thanks for the help of igor Rivin,m_goldberg, bbpdqn.

(Efficiency test)

IgorRivin = {{0, 1}, {1, 1}};
fibIgorRivin[n_] := MatrixPower[IgorRivin, n][[2, 1]]

fibIgorRivin[1000000] // AbsoluteTiming // First

0.2900166

fibGoldberg[n_Integer /; n > 0] := fibAux[1, 0, n]
fibAux[_, b_, 0] := b
fibAux[a_, b_, n_] := fibAux[a + b, a, n - 1]

 $IterationLimit::itlim: Iteration limit of 4096 exceeded. >>

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

fibUbpdqn[1000000] // AbsoluteTiming // First

 15.2288711

Fibonacci[1000000] // AbsoluteTiming // First

 0.2370136

Answer 1

A number of points:

First, if you need to compute several values of your sequence, your intial memo-ized implementation will NOT run into recursion limit problems.

Second, if you need to compute very few values, this method is extremely inefficient - the $n$-th Fibonacci number (or the $n$-th term of a linear recurrence in general) can be computed in $\log n$ operations. This can be done either by using linear algebra (the shift operator is an $n\times n$ matrix for an order $n$ recurrence) or the closed form solution (in terms of the roots of the characteristic polynomial.

For the matrix approach to Fibonacci:

fibmat = {{0, 1}, {1, 1}}
matfib[n_]:= MatrixPower[fibmat, n][[2, 1]]

matfib[200] // AbsoluteTiming
{0., 280571172992510140037611932413038677189525}

matfib[100000][[1]] // AbsoluteTiming
0.002000

As you see, this is much faster than the other ways if you need single values

Thirdly, once we are on the subject of closed form solutions, Mathematica does have RSolve[]/RSolveValue[], which work fine for linear recurrences in my experience.

Some words of explanation

First question is: how to transform a recurrence into matrix multiplication? well, if you have the recurrence (for simplicity I use order 2, but this obviously works for any order):

$$ a_n = c_1 a_{n-1} + c_2 a_{n-2}$$ then moving forward in "time" $n$ transforms the pair $a_{n-1}, a_n$ into $a_n, a_{n+1} = c_1 a_n + c_2 a_{n-1},$ so that this corresponds to acting on the vector $\begin{pmatrix} a_{n-1}\\a_n\end{pmatrix}$ by the matrix $$ \begin{pmatrix} 0 && 1\\ c_2 && c_1 \end{pmatrix}$$

So, computing the $n$th term of the recurrence corresponds to raising the matrix to the $n$th power (starting with the initial conditions).

Now, computing $x^n$ (where $x$ is a number or a matrix can be done by repeated squaring, and takes $O(\log n)$ multiplications. Of course, if you want to know the exact answer, the multiplications will take longer once you get further along, but this is not an issue when you are dealing with floating point numbers (you never increase the precision), and an amusing computation shows that this is a big win for exact integer computation IF asymptotically fast integer arithmetic is used (as it is in Mathematica).

I am 99% certain that this is the algorithm used for Fibonacci[] (of course it is even faster than what I did, since it is hardcoded in C, but if you have your own favorite recurrence, this is not going to be true.

Answer 2

If you can rewrite the recursion to be tail recursive, you will not run into recursion limits.

Here is an example of a tail-recursive implementation of factorial.

factorial[1, val_: 1] = val;
factorial[k_Integer /; k > 1, val_: 1] := factorial[k - 1, k val]

factorial /@ Range[10]

{1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800}

Block[{$RecursionLimit = 20}, factorial[42]]

1405006117752879898543142606244511569936384000000000

Note that in the last case, although the recursion limit was reduced to 20, no $RecursionLimit::reclim: message is emitted.

Update

By popular demand; i.e., belisarius' request.

fib[n_Integer /; n > 0] := fibAux[1, 0, n]
fibAux[_, b_, 0] := b
fibAux[a_, b_, n_] := fibAux[a + b, a, n - 1]

fib /@ Range[20]

{1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597, 2584, 4181, 6765}

Block[{$RecursionLimit = 20}, fib[200]]

280571172992510140037611932413038677189525

The general principle behind the two examples I've given is that the recursive function's argument sequence must carry all the state information needed at each step in the process.

---

References

http://en.wikipedia.org/wiki/Tail_call

What tools can help in realizing tail recursion?

Answer 3

Just approaches using Nest:

Fibonacci:

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

Note:

Timing[fi[200]]

yields:

{0., 280571172992510140037611932413038677189525}

or (really the same):

fb[1] := 1;
fb[2] := 1;
fb[n_] := First@Nest[{{1, 1}, {1, 0}}.# &, {1, 1}, n - 2]

yields same

Or factorial:

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

(noting for all these Mathematica has implementations or you could use e.g. Times@@Range[n] and define for n=0.

Answer 4

My thought about recursion formula this two weeks:

• For recursion formula $a_{n}=c_1a_{n-1}+c_2a_{n-2}\quad (n \geq 2)$

Thanks for @IgorRivin's explanation

$$ \begin{bmatrix} a_{n-1}\\ a_{n} \end{bmatrix} = \begin{bmatrix} a_{n-1}\\ c_1a_{n-1}+c_2a_{n-2} \end{bmatrix} = \begin{bmatrix} 0 & 1\\ c_2 & c_1 \end{bmatrix} \begin{bmatrix} a_{n-2}\\ a_{n-1} \end{bmatrix} $$ $$ =\ldots = \begin{bmatrix} 0 & 1\\ c_2 & c_1 \end{bmatrix}^{n-2} \begin{bmatrix} a_{1}\\ a_{2} \end{bmatrix} =A^{n-2}\beta $$

So for fib(n)=fib(n-1)+fib(n-2):

fib[n_Integer]:=
  Last@(MatrixPower[{{0, 1}, {1, 1}}, n - 2].{1, 1})

Test:

Fibonacci[100000]; // AbsoluteTiming

{0.0030002, Null}

fib[100000]; // AbsoluteTiming

{0.0040002, Null}

matfib[100000]; // AbsoluteTiming (*Igor Rivin*)

{0.0030002, Null}

• deCasteljau algorithm

  $$\vec{P}_{k,i}(u_0)=\color{red}{(1-u_0)}\vec{P}_{k-1,i}(u_0)+\color{red}{(u_0)}\vec{P}_{k-1,i+1}(u_0)$$

I discovered occationly a built-in function MovingAverage which can realize the deCasteljau algorithm

deCasteljau[pts : {{_, } ..}, u0] := Prepend[ NestList[ MovingAverage[#, {2, 3}] &, MovingAverage[pts, {1 - u0, u0}], Length@pts -2],pts]

 deCasteljau[pts : {{_, _} ..}, u0_] := 
   NestList[MovingAverage[#, {1 - u0, u0}] &, pts, Length@pts - 1]

test:

deCasteljau[{{0, 0}, {2, 4}, {4, 5}, {6, 0}}, 2/5]

  {{{0, 0}, {2, 4}, {4, 5}, {6, 0}}, {{4/5, 8/5}, {14/5, 22/5}, {24/5, 3}}, 
     {{8/5, 68/25}, {18/5, 96/25}}, {{12/5, 396/125}}}

This method can avoid using P[k_, i_, u0_] := (1 - u0) P[k - 1, i, u0] + u0 P[k - 1, i + 1, u0]

• Horner algorithm

Using Fold

 horner[coeffi_, x_Symbol] := 
   x #[[1]] + #[[2]] &@
     Fold[{x #[[1]] + #[[2]], #2} &, Take[coeffi, -2], Reverse@Drop[coeffi, -2]]

Test:

horner[{a0, a1, a2, a3, a4}, x]

 a0 + x (a1 + x (a2 + x (a4 + a3 x)))

Summarize

Sometimes, for me, using recursion formula in Mathematica can easily lead to Mathematica no responding(in this condition, I cannot use Alt+., command to abort caculation)

So it is necesary for me to try to using constructional operation to replace recursion formula.

Answer 5

Since the Fibonacci sequence is your example you should see:

• Fibonacci Sequence Generator

As described there you can also make use of LinearRecurrence:

LinearRecurrence[{1, 1}, {1, 1}, 15]

{1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610}

The syntax is simple and fairly general. It is comparatively fast for generating a complete sequence but it is not fast for finding single values.

In this earlier Stack Overflow Q&A Daniel Lichtblau gave the MatrixPower method that Igor Rivin describes. At that time I wrote a function to convert between the LinearRecurrence syntax and the MatrixPower equivalent:

linRec[r_, c_, i_] /; i <= Length@c := c[[i]]

linRec[r_, c_, i_] /; Length[r] == Length[c] :=
  With[{kern = Rest@IdentityMatrix[Length@r] ~Join~ {Reverse@r}},
    Last[MatrixPower[kern, i - Length@r]].c
  ]

Now:

linRec[{1, 1}, {1, 1}, 15]

610

And also:

LinearRecurrence[{-1, -1, 2, 3, 2}, {6, 8, 9, -9, -3}, {500}]

linRec[{-1, -1, 2, 3, 2}, {6, 8, 9, -9, -3}, 500]

{12519024244418274297611984359547397636366092718206074485347827784856486463673}

12519024244418274297611984359547397636366092718206074485347827784856486463673

The MatrixPower method is far faster than LinearRecurrence for a single value when the kernel is small:

LinearRecurrence[{-1, -1, 2, 3, 2}, {6, 8, 9, -9, -3}, {500000}] // AbsoluteTiming // First
linRec[{-1, -1, 2, 3, 2}, {6, 8, 9, -9, -3}, 500000]             // AbsoluteTiming // First

6.290360

0.030002

However you should be aware that with very large kernels it is slower than LinearRecurrence:

SeedRandom[1];
a = RandomInteger[{-9, 9}, {2, 200}];  (* kernel length 200 *)

LinearRecurrence[##, {500}] & @@ a // AbsoluteTiming // First
linRec[##, 500] & @@ a             // AbsoluteTiming // First

0.049003

0.362021

In such cases the Nest method I proposed in that old StackExchange Q&A is better:

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

lr[##, 500] & @@ a // AbsoluteTiming // First

0.017001

You can also generate the full sequence by using NestList if desired.

---

Related:

• Can the general term my recurrence equations be written with Floor or Mod?

• How can I get the general term of this recurrence equations?

---

Addendum

You asked about the difference between generating a complete sequence and finding single values.

• In some cases it is possible to directly calculate a single value, bypassing the recursion entirely. For example a Fibonacci number may be calculated using:

    Round[(GoldenRatio^(2 #) - (-1)^#)/(Sqrt[5] GoldenRatio^#)] &
  

  This is fast (though not as fast as MatrixPower) for a single value, but it would be very slow to calculate every value in the sequence using this function.

• The use of MatrixPower may appear iterative but I do not believe this is actually the case. Although I am not familiar with the algorithm used I would not expect all intermediate values to be generated in its calculation, therefore it is unsuited to generating the full sequence.

• Even if a method such as Nest ostensibly generates all the values in the sequence up to n the overhead (and memory consumption) of accumulating these values may larger than the computation itself, therefore it would be foolish to use e.g. NestList and then extract the last value.