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-1Mn = ( U Λ U-1) n
    = U Λ U-1U Λ U-1U Λ U-1 ... n times
    = U Λ Λ Λ ... U-1 
    = U Λ nU-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 = λ1X1X1,1 + X1,2 = λ1X1,1X1,1      = λ1X1,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.


The nth Fibonacci number is given by

f(n) = Floor(phi^n / sqrt(5) + 1/2) 

where

phi = (1 + sqrt(5)) / 2

Assuming that the primitive mathematical operations (+, -, * and /) are O(1) you can use this result to compute the nth Fibonacci number in O(log n) time (O(log n) because of the exponentiation in the formula).

In C#:

static double inverseSqrt5 = 1 / Math.Sqrt(5);
static double phi = (1 + Math.Sqrt(5)) / 2;
/* should use 
   const double inverseSqrt5 = 0.44721359549995793928183473374626
   const double phi = 1.6180339887498948482045868343656
*/

static int Fibonacci(int n) {
    return (int)Math.Floor(Math.Pow(phi, n) * inverseSqrt5 + 0.5);
}

If you want the exact number (which is a "bignum", rather than an int/float), then I'm afraid that

It's impossible!

As stated above, the formula for Fibonacci numbers is:

fib n = floor (phin/√5 + 1/2)

fib n ~= phin/√5

How many digits is fib n?

numDigits (fib n) = log (fib n) = log (phin/√5) = log phin - log √5 = n * log phi - log √5

numDigits (fib n) = n * const + const

it's O(n)

Since the requested result is of O(n), it can't be calculated in less than O(n) time.

If you only want the lower digits of the answer, then it is possible to calculate in sub-linear time using the matrix exponentiation method.


One of the exercises in SICP is about this, which has the answer described here.

In the imperative style, the program would look something like

Function Fib(count)
    a ← 1
    b ← 0
    p ← 0
    q ← 1

    While count > 0 Do
        If Even(count) Then
             pp² + q²
             q ← 2pq + q²
             countcount ÷ 2
        Else
             abq + aq + ap
             bbp + aq
             countcount - 1
        End If
    End While

    Return b
End Function