Following from Pillsy’s reference to matrix exponentiation, such that for the matrix
M = [1 1]
[1 0]
then
fib(n) = Mn1,2
Raising matrices to powers using repeated multiplication is not very efficient.
Two approaches to matrix exponentiation are divide and conquer which yields Mn in O(ln n) steps, or eigenvalue decomposition which is constant time, but may introduce errors due to limited floating point precision.
If you want an exact value greater than the precision of your floating point implementation, you have to use the O ( ln n ) approach based on this relation:
Mn = (Mn/2)2 if n even = M·Mn-1 if n is odd
The eigenvalue decomposition on M finds two matrices U and Λ such that Λ is diagonal and
M = U Λ U-1
Mn = ( U Λ U-1) n
= U Λ U-1 U Λ U-1 U Λ U-1 ... n times
= U Λ Λ Λ ... U-1
= U Λ n U-1
Raising a the diagonal matrix Λ to the nth power is a simple matter of raising each element in Λ to the nth, so this gives an O(1) method of raising M to the nth power. However, the values in Λ are not likely to be integers, so some error will occur.
Defining Λ for our 2×2 matrix as
Λ = [ λ1 0 ] = [ 0 λ2 ]
To find each λ, we solve
|M - λI| = 0
which gives
|M - λI| = -λ ( 1 - λ ) - 1 λ² - λ - 1 = 0
using the quadratic formula
λ = ( -b ± √ ( b² - 4ac ) ) / 2a
= ( 1 ± √5 ) / 2
{ λ1, λ2 } = { Φ, 1-Φ } where Φ = ( 1 + √5 ) / 2
If you’ve read Jason’s answer, you can see where this is going to go.
Solving for the eigenvectors X1 and X2:
if X1 = [ X1,1, X1,2 ] M.X1 1 = λ1X1 X1,1 + X1,2 = λ1 X1,1 X1,1 = λ1 X1,2 => X1 = [ Φ, 1 ] X2 = [ 1-Φ, 1 ]
These vectors give U:
U = [ X1,1, X2,2 ]
[ X1,1, X2,2 ]
= [ Φ, 1-Φ ]
[ 1, 1 ]
Inverting U using
A = [ a b ]
[ c d ]
=>
A-1 = ( 1 / |A| ) [ d -b ]
[ -c a ]
so U-1 is given by
U-1 = ( 1 / ( Φ - ( 1 - Φ ) ) [ 1 Φ-1 ]
[ -1 Φ ]
U-1 = ( √5 )-1 [ 1 Φ-1 ]
[ -1 Φ ]
Sanity check:
UΛU-1 = ( √5 )-1 [ Φ 1-Φ ] . [ Φ 0 ] . [ 1 Φ-1 ]
[ 1 1 ] [ 0 1-Φ ] [ -1 Φ ]
let Ψ = 1-Φ, the other eigenvalue
as Φ is a root of λ²-λ-1=0
so -ΨΦ = Φ²-Φ = 1
and Ψ+Φ = 1
UΛU-1 = ( √5 )-1 [ Φ Ψ ] . [ Φ 0 ] . [ 1 -Ψ ]
[ 1 1 ] [ 0 Ψ ] [ -1 Φ ]
= ( √5 )-1 [ Φ Ψ ] . [ Φ -ΨΦ ]
[ 1 1 ] [ -Ψ ΨΦ ]
= ( √5 )-1 [ Φ Ψ ] . [ Φ 1 ]
[ 1 1 ] [ -Ψ -1 ]
= ( √5 )-1 [ Φ²-Ψ² Φ-Ψ ]
[ Φ-Ψ 0 ]
= [ Φ+Ψ 1 ]
[ 1 0 ]
= [ 1 1 ]
[ 1 0 ]
= M
So the sanity check holds.
Now we have everything we need to calculate Mn1,2:
Mn = UΛnU-1
= ( √5 )-1 [ Φ Ψ ] . [ Φn 0 ] . [ 1 -Ψ ]
[ 1 1 ] [ 0 Ψn ] [ -1 Φ ]
= ( √5 )-1 [ Φ Ψ ] . [ Φn -ΨΦn ]
[ 1 1 ] [ -Ψn ΨnΦ ]
= ( √5 )-1 [ Φ Ψ ] . [ Φn Φn-1 ]
[ 1 1 ] [ -Ψn -Ψn-1 ] as ΨΦ = -1
= ( √5 )-1 [ Φn+1-Ψn+1 Φn-Ψn ]
[ Φn-Ψn Φn-1-Ψn-1 ]
so
fib(n) = Mn1,2
= ( Φn - (1-Φ)n ) / √5
Which agrees with the formula given elsewhere.
You can derive it from a recurrance relation, but in engineering computing and simulation calculating the eigenvalues and eigenvectors of large matrices is an important activity, as it gives stability and harmonics of systems of equations, as well as allowing raising matrices to high powers efficiently.