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