# 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.
#
