nth fibonacci number in sublinear time

Biswajyoti Das picture Biswajyoti Das · Oct 6, 2009 · Viewed 39.3k times · Source

Is there any algorithm to compute the nth fibonacci number in sub linear time?

Answer

Pete Kirkham picture Pete Kirkham · Oct 6, 2009

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 2x2 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 ]   [ -Ψnn-1 ] as ΨΦ = -1

   = ( √5 )-1 [ Φn+1n+1      Φnn ]
              [ Φnn      Φn-1n-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.