# 1.30 Lesson 30

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