How to deal with recursion formula in Mathematica?

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
• You could simply set the RecursionLimit to a higher value. For example $RecursionLimit = 1000 – paw Oct 1 '14 at 2:37 • With enough beer, you can tweak any recursive or iterative construct to do what you want. But why? The most reasonable structure to represent the Fib. progression is the recursive one. Oct 1 '14 at 2:42 • What version of Mathematica are you using? In V10 on OS X and, I think V9 as well, the default value for$RecursionLimit is 1024. Oct 1 '14 at 3:15
• I recommend that you look at the function RSolve which can determine closed form representations for many recursive equations. Oct 1 '14 at 5:58
• Have a look at this approach if it can help How to find the sum all even numbers of this sequence?. Oct 1 '14 at 7:40

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.

• +1,As the n becoming larger and larger ,like 1000000,your alogrithm is fastest.:-),BTW, could you give a explanation to your algorithm? Thanks !!
– xyz
Oct 1 '14 at 13:34
• @Tangshutao see the edit... Oct 2 '14 at 12:20

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? • Perhaps posting the tail recursive code for the Fibonacci numbers would be a beau geste towards the OP :).... (and me) Oct 1 '14 at 4:20 • +1, Powerful, fib[800] can be caculated immediately when$RecursionLimit = 256. I think I can apply this idea to other problem.:-) Thanks!
– xyz
Oct 1 '14 at 11:01

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.

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.

Since the Fibonacci sequence is your example you should see:

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:

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.

• I have a question: what's the difference between generating a complete sequence and finding single values?I'm always thinking that: 1,generating a complete sequence 2, Exacting the single values from the complete sequence.@Mr.Wizard, Could you explain it to me?Thanks sincerely!!:-)
– xyz
Oct 12 '14 at 6:06
• @Tangshutao I am too sleepy tonight but I will try to add a few comments on that tomorrow. Oct 12 '14 at 6:09
• I'm so sorry,@Mr.Wizard I forgoit the problem timezone beteen USA and China.
– xyz
Oct 12 '14 at 6:15
• @ShutaoTang Please see the Addendum. I hope that I understood your question and that you find my comments useful. Oct 15 '14 at 18:50