# 1.29 Lesson 29

```# Lesson 29.  Fibonacci Numbers
# -----------------------------
#
# One way to define a sequence a[n] is by a recurrence relation, which is an equation expressing a[n] in terms of
# the a[k] for k < n.  For example,
> a[n] = 2*a[n-1]:
# This alone doesn't determine the a[n]'s.  You also need an initial condition: in this case a value for a[0].
> a[0]:= 1;
> for nn from 1 to 3 do a[nn]:= 2*a[nn-1] od;

a[0] := 1

a[1] := 2

a[2] := 4

a[3] := 8

# Obviously in this case we'll have a[n] = 2^n for all n.  Here's a slightly more complicated relation:
> F[n] = F[n-1]+F[n-2]:
> F[0] := 0: F[1] := 1:
# This sequence is known as the Fibonacci numbers.  Note that here we need two initial values to determine the
# sequence.
> for nn from 2 to 10 do
>   F[nn]:= F[nn-1] + F[nn-2] od:
> seq(F[nn], nn=0..10);

0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55

# Here's a combinatorial problem where the Fibonacci numbers come up.
# Suppose there are n Math 210 classes in the term.
# You might skip some of the classes, but you don't want to ever miss two consecutive classes.  How many
# different attendance patterns can you have?  e.g. with n = 3, here are the 5 possible patterns (0 = skip, 1 =
# attend):
# 010, 011, 101, 110, 111
# You can get a recurrence relation this way.  Let A(n) be the number of patterns for n classes.   If you attend
# the first class, you have A(n-1) possible patterns for the remaining n-1 classes.  On the other hand, if you skip
# the first class, you must attend the second; after that you have A(n-2) possible patterns for the remaining n-2
# classes.
# So A(n) = A(n-1) + A(n-2).
# The initial conditions are A(1)=2 and A(0)=1 (or if you don't think A(0) should be defined, A(2)=3).  Note that
# A(n)=F[n+2] for n=0, 1, 2.  By  mathematical induction, A(n) = F[n+2] for all n.
#
# How can we calculate the Fibonacci numbers?  The most straightforward way is with a "for" loop as we had
# above.  To calculate F[n], you need to calculate F[k] for all k < n.
#
# Here's another way that looks neater.  Instead of explicitly calculating F[k] for 1 <= k < n, get Maple to do it
# automatically when needed.
> fib1:= n -> fib1(n-1) + fib1(n-2):
> fib1(0):= 0: fib1(1):= 1:
# This is taking advantage of the "remember table" to supply the values of fib1(0) and fib1(1) rather than using
# the recurrence for them.  Don't forget to define those values or you'll get an infinite recursion.  Also don't try
# fib1(n) where n is not a positive integer.
# But this is a really terrible way to calculate the Fibonacci numbers.
# Let T(n) be the number of steps "fib1(n)" takes to do the calculation, where each step is either an addition
# (according to the recursive formula) or remembering the value of fib1(0) or fib1(1).
# Then for n >= 2, T(n) = 1 + T(n-1) + T(n-2), with T(0) = 1 and T(1) = 1.  We can write the recurrence
# relation as
> T(n) + 1 = (T(n-1)+1) + (T(n-2)+1):
> T(0)+1 = 2: T(1)+1 = 2:
# From this it is easy to see that T(n) = 2 F[n+1] -1.
#
# A better way would be to have Maple remember the values of each Fibonacci number it calculates.  A Maple
# procedure will do this if you specify "option remember" in its definition.  You can't use "->" to define such a
# procedure: instead you use "proc":
> fib2:= proc(n) option remember;
>  fib2(n-1)+fib2(n-2);
> end;

fib2 := proc(n)
options remember;
fib2(n-1)+fib2(n-2)
end

> fib2(0):= 0: fib2(1):= 1:
> fib2(30);

832040

# To calculate fib2(n), Maple needs to calculate fib2(n-1), but then it just remembers fib2(n-2) instead of
# calculating it again.  So the number of steps is approximately 2 n.
# There's one problem you might encounter with "fib2" on your home computer (it won't be so bad in the lab).
> fib2(200);
Error, (in fib2) STACK OVERFLOW

# What does Maple do when you press "Enter" here?  To calculate fib2(200), according to the procedure
# definition, it must first find fib2(199), then fib2(198), then add them.  While it is calculating fib2(199), it must
# keep in the computer's memory some information that will let it complete the calculation of fib2(200) when
# fib2(199) finishes.  The area of memory where this is kept is called the stack, and the information is called a
# stack frame.  By the time it gets down to fib2(30), which is the highest value in the remember table, the stack
# has frames for fib2(200), fib2(199), ..., fib2(31).  The stack size is rather limited in PC and Mac versions of
# Maple, and this leads to the stack overflow message.
#
# The next methods of calculating Fibonacci numbers require some linear algebra.  We can express the
# recurrence relation in terms of vectors and matrices.  Let X(n) be the vector [F[n-1], F[n]], so X(n+1) = [F[n],
# F[n+1]] = [F[n], F[n-1] + F[n]].  This can be
# written as X(n+1) = M X(n) for a certain matrix M.
> with(linalg): M:= matrix([[0,1],[1,1]]);
Warning: new definition for   norm
Warning: new definition for   trace

[ 0  1 ]
M := [      ]
[ 1  1 ]

> [F[n],F[n+1]] = evalm(M &* [F[n-1],F[n]]);

[F[n], F[n + 1]] = [ F[n], F[n - 1] + F[n] ]

# Starting from X(1) = [F[0], F[1]] = [0,1], we have X(n) = M^(n-1) X(1).  F[n] is the 2nd component of X(n),
# which is the [2,2] matrix element of M^(n-1).
> fib3:= proc(n) option remember;
>     evalm(M^(n-1))[2,2] end:
> fib3(200);

280571172992510140037611932413038677189525

# This procedure hides the details of the calculation in the workings of "evalm".  What's an efficient way to
# calculate the n'th power of a matrix (other than by multiplying that many copies of the matrix together)?  The
# method is called "repeated squaring".  Thus to calculate M^64, we would first calculate M^2, then square that
# to get M^4, then square that to get M^8, then M^16, M^32, and M^64.  It only required 6 squarings.  To
# calculate M^199, we write 199 as a sum of powers of 2 (this is the base-2 representation of 199):
#  199 = 1 + 2 + 4 + 64 + 128
# so M^199 = M  M^2 M^4 M^64 M^128.
# A convenient way to program it is the following:
> mpow:= proc(n) option remember;\
if type(n,even) then \
evalm(mpow(n/2)^2)\
elif type(n,odd) then\
evalm(mpow((n-1)/2)^2 &* M)\
else M^n \
fi\
end:
> mpow(0):= matrix([[1,0],[0,1]]):
> mpow(1):= M:
> mpow(200);

[173402521172797813159685037284371942044301,

280571172992510140037611932413038677189525]

[280571172992510140037611932413038677189525,

453973694165307953197296969697410619233826]

> fib4:= n -> mpow(n-1)[2,2];

fib4 := n -> mpow(n - 1)[2, 2]

> fib4(1000);

434665576869374564356885276750406258025646605173717804024817290\
895365554179490518904038798400792551692959225930803226347752096\
896232398733224711616429964409065331879382989696499285160037044\
76137795166849228875

# As it stands, "fib4" works only for positive integers.  We'd like something that works also for symbolic
# expressions, e.g we'd like to be able to get F[2*n+3]  as an expression in terms of F[n] and F[n-1].  Note that
# M^n is the matrix
# [ F[n-1]     F[n]   ]
# [ F[n]       F[n+1] ].
> evalm(M^4);

[ 2  3 ]
[      ]
[ 3  5 ]

> mpow:= proc(n) option remember;
>   if type(n,posint) then
>     if type(n,even) then evalm(mpow(n/2)^2)
>     else evalm(mpow((n-1)/2)^2 &* M)
>     fi
>   elif type(n,negint) then
>     evalm(mpow(-n)^(-1))
>   elif type(n,`+`) then
>     evalm(`&*`(op(map(mpow,n))))
>   elif type(n, `*`) then
>     evalm(mpow(n/op(1,n))^op(1,n))
>   else matrix([[F[n-1],F[n]], [F[n],F[n]+F[n-1]]])
>   fi
> end:
> mpow(0):= matrix([[1,0],[0,1]]):
> fib4(n);

F[n]

> F[n+3] = fib4(n+3);

F[n + 3] = 3 F[n] + 2 F[n - 1]

> F[n-2] = fib4(n-2);

F[n - 2] = - F[n - 1] + F[n]

> F[n+m] = fib4(n+m);

F[n + m] = F[n - 1] F[m] + F[n] F[m] + F[n] F[m - 1]

# Where did this come from?
> M^(m+n) = mpow(m) &* mpow(n);

(m + n)
M        =

[ F[m - 1]        F[m]      ]    [ F[n - 1]        F[n]      ]
[                           ] &* [                           ]
[   F[m]    F[m] + F[m - 1] ]    [   F[n]    F[n] + F[n - 1] ]

# This was only one of many interesting identities involving the Fibonacci numbers.
# There is a whole journal, the Fibonacci Quarterly, devoted to Fibonacci numbers
# and related recurrence relations.  Near the back of each issue is a problem section, containing problems at
# various levels of difficulty contributed by readers.  Many of
# the problems in the elementary section can be solved quite readily using "fib4".
# Here are some examples:
#
# (1) Let H[n] be any solution of the recurrence H[n] = H[n-1] + H[n-2].  Show
# that 7 H[n] = H[n+15] mod 10.
# ("=" here should be "is congruent to", which is written with three horizontal lines,
# but I don't think I can write it in this font.  In any case, this means that 7 H[n] and
# H[n+15] have the same remainder when divided by 10, i.e. the same last digit in
# the decimal representation.
#
# Solution: If X(n) = [ H[n], H[n+1]], then X(n+1) = M X(n), so X(n+15) = M^15 X(n).
> X:= [ H[n], H[n+1]]:
> evalm(M^15 &* X - 7*X);

[ 370 H[n] + 610 H[n + 1], 610 H[n] + 980 H[n + 1] ]

# The first component is H[n+15] - 7 H[n], and clearly it is divisible by 10.
#
# (2) Show that F[5^n] is divisible by 5^n.
#
# Solution:
> fib4(5*k);

4              2     3         5
5 F[k] F[k - 1]  + 20 F[k - 1]  F[k]  + 5 F[k]

4              3     2
+ 15 F[k - 1] F[k]  + 10 F[k - 1]  F[k]

> "/(5*F[k]);

4              2     3         5
1/5 (5 F[k] F[k - 1]  + 20 F[k - 1]  F[k]  + 5 F[k]

4              3     2
+ 15 F[k - 1] F[k]  + 10 F[k - 1]  F[k] )/F[k]

> expand(");

4             2     2       4                  3
F[k - 1]  + 4 F[k - 1]  F[k]  + F[k]  + 3 F[k - 1] F[k]

3
+ 2 F[k - 1]  F[k]

# Since F[5*k] is divisible by 5 F[k], we get that F[5^n] is divisible by 5^n for every n.
#

```