# Lesson 30. Fibonacci Numbers, continued
# ----------------------------------------
#
# Last time we looked at some ways of calculating the Fibonacci numbers. I'll just use the last one, which I'll call
# "fib" instead of "fib4".
> restart: with(linalg):
Warning: new definition for norm
Warning: new definition for trace
> M:= matrix([[0,1],[1,1]]);
[ 0 1 ]
M := [ ]
[ 1 1 ]
> 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:
> fib:= n -> mpow(n-1)[2,2];
fib := n -> mpow(n - 1)[2, 2]
> fib(1000);
4346655768693745643568852767504062580256466051737178040248172908953655541794905\
1890403879840079255169295922593080322634775209689623239873322471161642996440906\
533187938298969649928516003704476137795166849228875
# As it stands, "fib" 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]]):
# I have saved these definitions (and the "with(linalg)" command) in a file "fib.def" and placed it in the "doc"
# directory. You can read them into your own worksheets with the command
> read `doc/fib.def`;
> fib(n);
F[n]
> F[n+3] = fib(n+3);
F[n + 3] = 3 F[n] + 2 F[n - 1]
> F[n-2] = fib(n-2);
F[n - 2] = - F[n - 1] + F[n]
> F[n+m] = fib(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);
(n + m) [ F[m - 1] F[m] ] [ F[n - 1] F[n] ]
M = [ ] &* [ ]
[ 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 "mpow" and "fib".
# 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 can't 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 (but not by 5^(n+1)).
#
# Solution (to first part - I won't bother with the second):
> fib(5*k);
4 2 3 5 4
5 F[k] F[k - 1] + 20 F[k - 1] F[k] + 5 F[k] + 15 F[k - 1] F[k]
3 2
+ 10 F[k - 1] F[k]
> "/(5*F[k]);
4 2 3 5 4
1/5 (5 F[k] F[k - 1] + 20 F[k - 1] F[k] + 5 F[k] + 15 F[k - 1] F[k]
3 2
+ 10 F[k - 1] F[k] )/F[k]
> expand(");
4 2 2 4 3 3
F[k - 1] + 4 F[k - 1] F[k] + F[k] + 3 F[k - 1] F[k] + 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. This is mathematical
# induction again:
> F[5^0] = fib(5^0);
F[1] = 1
# k = 0: F[5^0] is divisible by 5^0.
# Induction step: if F[5^k] is divisible by 5^k, then F[5^(k+1)] = F[5*5^k] is divisible by 5 F[5^k], and thus by
# 5^(k+1).
#
# (3) What is a[n], if a[n+1] = 5 a[n]^3 - 3 a[n] for n >= 0 with a[0] = 1?
#
# This is another recurrence relation, but a nonlinear one. Let's start
# by calculating a few terms.
> a[0]:= 1;
a[0] := 1
> for nn from 0 to 5 do \
a[nn+1]:= 5*a[nn]^3 - 3*a[nn] od;
a[1] := 2
a[2] := 34
a[3] := 196418
a[4] := 37889062373143906
a[5] := 271964099255182923543922814194423915162591622175362
a[6] :=
1005784040477634373949250581829128589677588166452408546692608465300231471006495\
79758906802170865134804283198088355735627634798716539277515304438631163554
# Well, it may be good that we stopped here: the numbers are growing very rapidly. This being a problem that's
# supposed to have something to do with Fibonacci numbers, we might remember that a[2]=34 is a Fibonacci
# number:
> fib(9);
34
# Of course a[1] = 2 is also a Fibonacci number, F[3], and a[1] = 1 is F[2] and F[1]. This might lead you to make
# a wild guess that a[n] = F[3^n]. Try it for n=3,4,5, and6:
> seq(fib(3^kk)-a[kk], kk = 3 .. 6);
0, 0, 0, 0
# It seems to work. Now can we prove it?
> fib(3*k);
2 3 2
3 F[k] F[k - 1] + 2 F[k] + 3 F[k - 1] F[k]
# This should be 5 F[k]^3 - 3 F[k], at least if k is a power of 3, in order for our conjecture to be true.
> difference:= " - (5*F[k]^3 - 3*F[k]);
2 3 2
difference := 3 F[k] F[k - 1] - 3 F[k] + 3 F[k - 1] F[k] + 3 F[k]
> factor(");
2 2
3 F[k] (F[k - 1] - F[k] + F[k - 1] F[k] + 1)
# Well, the 3 F[k] isn't 0, so maybe the other factor is.
> Q:= k -> fib(k-1)^2 - fib(k)^2 + fib(k-1)*fib(k)+1;
2 2
Q := k -> fib(k - 1) - fib(k) + fib(k - 1) fib(k) + 1
> seq(Q(kk), kk = 0 .. 5);
2, 0, 2, 0, 2, 0
# It looks like this should be 0 for odd k. That will be fine for our purposes: the powers of 3 are all odd. Can we
# prove it by mathematical induction?
> Q(k+1);
2 2
F[k] - (F[k] + F[k - 1]) + F[k] (F[k] + F[k - 1]) + 1
# If the pattern works, Q(k+1) should be 2 - Q(k), so 2 for even k becomes 0 for odd k and vice versa.
> expand(") - (2 - Q(k));
0
# Yes, it works: Q(k+1) = 2 - Q(k).
# There's another way to look at this result about Q, which can be rewritten as
# F[k-1]^2 - F[k]^2 + F[k-1] F[k] = (-1)^k. The left side is the determinant of the matrix
# M^k = [ F[k-1] F[k] ] = [ F[k-1] F[k] ]
# [ F[k] F[k+1] ] [ F[k] F[k-1]+F[k] ].
# But det(M^k) = (det M)^k = (-1)^k.
#
