nth fibonacci number in sublinear time

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

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 Mⁿ 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 X₁ and X₂:

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   Φ ]
   = ( &radic;5 )<sup>-1</sup> [ &Phi;   &Psi; ] . [ &Phi;   -&Psi;&Phi; ] 
             [ 1   1 ]   [ -&Psi;  &Psi;&Phi; ]

   = ( &radic;5 )<sup>-1</sup> [ &Phi;   &Psi; ] . [ &Phi;    1 ] 
             [ 1   1 ]   [ -&Psi;  -1 ]

   = ( &radic;5 )<sup>-1</sup> [ &Phi;&sup2;-&Psi;&sup2;  &Phi;-&Psi; ] 
              [ &Phi;-&Psi;      0 ]

   = [ &Phi;+&Psi;   1 ]    
     [ 1     0 ]

   = [ 1     1 ] 
     [ 1     0 ]

   = <b>M</b> 

So the sanity check holds.

Now we have everything we need to calculate Mⁿ_1,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.

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 (phiⁿ/√5 + ¹/₂) > > fib n ~= phiⁿ/√5

How many digits is fib n?

> numDigits (fib n) = log (fib n) = log (phiⁿ/√5) = log phiⁿ - 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 http://mitpress.mit.edu/sicp/full-text/book/book-Z-H-11.html#%_thm_1.19">exercises in SICP is about this, which has the answer described http://www.kendyck.com/archives/2005/05/13/solution-to-sicp-exercise-119/">here.</a>

In the imperative style, the program would look something like

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

<b>While</b> <i>count</i> &gt; 0 <b>Do</b>
    <b>If</b> Even(<i>count</i>) <b>Then</b>
         <i>p</i> ← <i>p</i>² + <i>q</i>²
         <i>q</i> ← 2<i>pq</i> + <i>q</i>²
         <i>count</i> ← <i>count</i> ÷ 2
    <b>Else</b>
         <i>a</i> ← <i>bq</i> + <i>aq</i> + <i>ap</i>
         <i>b</i> ← <i>bp</i> + <i>aq</i>
         <i>count</i> ← <i>count</i> - 1
    <b>End If</b>
<b>End While</b>

<b>Return</b> <i>b</i>

End Function

You can do it by exponentiating a matrix of integers as well. If you have the matrix

    / 1  1 \
M = |      |
    \ 1  0 /

then (M^n)[1, 2] is going to be equal to the nth Fibonacci number, if [] is a matrix subscript and ^ is matrix exponentiation. For a fixed-size matrix, exponentiation to an positive integral power can be done in O(log n) time in the same way as with real numbers.

EDIT: Of course, depending on the type of answer you want, you may be able to get away with a constant-time algorithm. Like the other formulas show, the nth Fibonacci number grows exponentially with n. Even with 64-bit unsigned integers, you'll only need a 94-entry lookup table in order to cover the entire range.

SECOND EDIT: Doing the matrix exponential with an eigendecomposition first is exactly equivalent to JDunkerly's solution below. The eigenvalues of this matrix are the (1 + sqrt(5))/2 and (1 - sqrt(5))/2.

Wikipedia has a closed form solution http://en.wikipedia.org/wiki/Fibonacci_number

Or in c#:

    public static int Fibonacci(int N)
    {
        double sqrt5 = Math.Sqrt(5);
        double phi = (1 + sqrt5) / 2.0;
        double fn = (Math.Pow(phi, N) - Math.Pow(1 - phi, N)) / sqrt5;
        return (int)fn;
    }

For really big ones, this recursive function works. It uses the following equations:

F(2n-1) = F(n-1)^2 + F(n)^2
F(2n) = (2*F(n-1) + F(n)) * F(n)

You need a library that lets you work with big integers. I use the BigInteger library from https://mattmccutchen.net/bigint/.

Start with an array of of fibonacci numbers. Use fibs[0]=0, fibs[1]=1, fibs[2]=1, fibs[3]=2, fibs[4]=3, etc. In this example, I use an array of the first 501 (counting 0). You can find the first 500 non-zero Fibonacci numbers here: http://home.hiwaay.net/~jalison/Fib500.html. It takes a little editing to put it in the right format, but that is not too hard.

Then you can find any Fibonacci number using this function (in C):

BigUnsigned GetFib(int numfib)
{
int n;
BigUnsigned x, y, fib;  

if (numfib < 501) // Just get the Fibonacci number from the fibs array
    {
       fib=(stringToBigUnsigned(fibs[numfib]));
    }
else if (numfib%2) // numfib is odd
    {
       n=(numfib+1)/2;
       x=GetFib(n-1);
       y=GetFib(n);
       fib=((x*x)+(y*y));
    }
else // numfib is even
    {
       n=numfib/2;
       x=GetFib(n-1);
       y=GetFib(n);
       fib=(((big2*x)+y)*y);
   }
return(fib);
}

I've tested this for the 25,000th Fibonacci number and the like.

Here's my recursive version that recurses log(n) times. I think that it's easiest to read in the recursive form:

def my_fib(x):
  if x < 2:
    return x
  else:
    return my_fib_helper(x)[0]

def my_fib_helper(x):
  if x == 1:
    return (1, 0)
  if x % 2 == 1:
    (p,q) = my_fib_helper(x-1)
    return (p+q,p)
  else:
    (p,q) = my_fib_helper(x/2)
    return (p*p+2*p*q,p*p+q*q)

It works because you can compute fib(n),fib(n-1) using fib(n-1),fib(n-2) if n is odd and if n is even, you can compute fib(n),fib(n-1) using fib(n/2),fib(n/2-1).

The base case and the odd case are simple. To derive the even case, start with a,b,c as consecutive fibonacci values (eg, 8,5,3) and write them in a matrix, with a = b+c. Notice:

[1 1] * [a b]  =  [a+b a]
[1 0]   [b c]     [a   b]

From that, we see that a matrix of the first three fibonacci numbers, times a matrix of any three consecutive fibonacci numbers, equals the next. So we know that:

      n
[1 1]   =  [fib(n+1) fib(n)  ]
[1 0]      [fib(n)   fib(n-1)]

So:

      2n                        2
[1 1]    =  [fib(n+1) fib(n)  ]
[1 0]       [fib(n)   fib(n-1)]

Simplifying the right hand side leads to the even case.

using R

l1 <- (1+sqrt(5))/2
l2 <- (1-sqrt(5))/2

P <- matrix(c(0,1,1,0),nrow=2) #permutation matrix
S <- matrix(c(l1,1,l2,1),nrow=2)
L <- matrix(c(l1,0,0,l2),nrow=2)
C <- c(-1/(l2-l1),1/(l2-l1))

k<-20 ; (S %*% L^k %*% C)[2]
[1] 6765

Fixed point arithmetic is inaccurate. Jason's C# code gives incorrect answer for n = 71 (308061521170130 instead of 308061521170129) and beyond.

For correct answer, use a computational algebra system. Sympy is such a library for Python. There's an interactive console at http://live.sympy.org/ . Copy and paste this function

phi = (1 + sqrt(5)) / 2
def f(n):
    return floor(phi**n / sqrt(5) + 1/2)

Then calculate

>>> f(10)
55

>>> f(71)
308061521170129

You might like to try inspecting phi.

Here's a one-liner that computes F(n), using integers of size O(n), in O(log n) arithmetic operations:

for i in range(1, 50):
	print(i, pow(2<<i, i, (4<<2*i)-(2<<i)-1)//(2<<i))

Using integers of size O(n) is reasonable, since that's comparable to size of the answer.

To understand this, let phi be the golden ratio (the largest solution to x^2=x+1) and F(n) be the n'th Fibonacci number, where F(0)=0, F(1)=F(2)=1

Now, phi^n = F(n-1) + F(n)phi.

> Proof by induction: phi^1 = 0 + 1*phi = F(0) + F(1)phi. And if phi^n = > F(n-1) + F(n)phi, then phi^(n+1) = F(n-1)phi + F(n)phi^2 = F(n-1)phi + > F(n)(phi+1) = F(n) + (F(n)+F(n-1))phi = F(n) + F(n+1)phi. The only tricky step in this calculation is the one that replaces phi^2 by (1+phi), which follows because phi is the golden ratio.

Also numbers of the form (a+b*phi), where a, b are integers are closed under multiplication.

> Proof: (p0+p1phi)(q0+q1phi) = p0q0 + (p0q1+q1p0)phi + p1q1phi^2 = > p0q0 + (p0q1+q1p0)phi + p1q1(phi+1) = (p0q0+p1q1) + > (p0q1+q1p0+p1q1)*phi.

Using this representation, one can compute phi^n in O(log n) integer operations using exponentiation by squaring. The result will be F(n-1)+F(n)phi, from which one can read off the n'th Fibonacci number.

def mul(p, q):
	return p[0]*q[0]+p[1]*q[1], p[0]*q[1]+p[1]*q[0]+p[1]*q[1]

def pow(p, n):
	r=1,0
	while n:
		if n&1: r=mul(r, p)
		p=mul(p, p)
		n=n>>1
	return r

for i in range(1, 50):
	print(i, pow((0, 1), i)[1])

Note that the majority of this code is a standard exponentiation-by-squaring function.

To get to the one-liner that starts this answer, one can note that representing phi by a large enough integer X, one can perform (a+b*phi)(c+d*phi) as the integer operation (a+bX)(c+dX) modulo (X^2-X-1). Then the pow function can be replaced by the standard Python pow function (which conveniently includes a third argument z which calculates the result modulo z. The X chosen is 2<<i.

Apart from fine-tuning by mathematical approaches, one of the best optimum solution (I believe) is using a dictionary in order to avoid repetitive calculations.

import time

_dict = {1:1, 2:1}

def F(n, _dict):
    if n in _dict.keys():
        return _dict[n]
    else:
        result = F(n-1, _dict) + F(n-2, _dict)
        _dict.update({n:result})
        return result

start = time.time()

for n in range(1,100000):
    result = F(n, _dict) 

finish = time.time()

print(str(finish - start))

We start with trivial dictionary (first two values of Fibonacci sequence) and constantly adding Fibonacci values to the dictionary.

It took about 0.7 seconds for the first 100000 Fibonacci values (Intel Xeon CPU E5-2680 @ 2.70 GHz, 16 GB RAM, Windows 10-64 bit OS)

see divide and conquer algorithm here

The link has pseudocode for the matrix exponentiation mentioned in some of the other answers for this question.

You can use the weird square rooty equation to get an exact answer. The reason is that the $\sqrt(5)$ falls out at the end, you just have to keep track of the coefficients with your own multiplication format.

def rootiply(a1,b1,a2,b2,c):
    ''' multipy a1+b1*sqrt(c) and a2+b2*sqrt(c)... return a,b'''
    return a1*a2 + b1*b2*c, a1*b2 + a2*b1

def rootipower(a,b,c,n):
    ''' raise a + b * sqrt(c) to the nth power... returns the new a,b and c of the result in the same format'''
    ar,br = 1,0
    while n != 0:
        if n%2:
            ar,br = rootiply(ar,br,a,b,c)
        a,b = rootiply(a,b,a,b,c)
        n /= 2
    return ar,br

def fib(k):
    ''' the kth fibonacci number'''
    a1,b1 = rootipower(1,1,5,k)
    a2,b2 = rootipower(1,-1,5,k)
    a = a1-a2
    b = b1-b2
    a,b = rootiply(0,1,a,b,5)
    # b should be 0!
    assert b == 0
    return a/2**k/5

if __name__ == "__main__":
    assert rootipower(1,2,3,3) == (37,30) # 1+2sqrt(3) **3 => 13 + 4sqrt(3) => 39 + 30sqrt(3)
    assert fib(10)==55

I have come across some of the methods for calculating Fibonacci with efficient time complexity following are some of them -

Method 1 - Dynamic Programming Now here the substructure is commonly known hence I'll straightly Jump to the solution -

static int fib(int n) 
{ 
int f[] = new int[n+2]; // 1 extra to handle case, n = 0 
int i; 

f[0] = 0; 
f[1] = 1; 
  
for (i = 2; i <= n; i++) 
{ 
    f[i] = f[i-1] + f[i-2]; 
} 
   
return f[n]; 
}

A space-optimized version of above can be done as follows -

static int fib(int n) 
 { 
    int a = 0, b = 1, c; 
    if (n == 0) 
        return a; 
    for (int i = 2; i <= n; i++) 
    { 
        c = a + b; 
        a = b; 
        b = c; 
    } 
    return b; 
} 

Method 2- ( Using power of the matrix {{1,1},{1,0}} )

This an O(n) which relies on the fact that if we n times multiply the matrix M = {{1,1},{1,0}} to itself (in other words calculate power(M, n )), then we get the (n+1)th Fibonacci number as the element at row and column (0, 0) in the resultant matrix. This solution would have O(n) time.

The matrix representation gives the following closed expression for the Fibonacci numbers: fibonaccimatrix

static int fib(int n) 
{ 
int F[][] = new int[][]{{1,1},{1,0}}; 
if (n == 0) 
	return 0; 
power(F, n-1); 

return F[0][0]; 
} 

/*multiplies 2 matrices F and M of size 2*2, and 
puts the multiplication result back to F[][] */
static void multiply(int F[][], int M[][]) 
{ 
int x = F[0][0]*M[0][0] + F[0][1]*M[1][0]; 
int y = F[0][0]*M[0][1] + F[0][1]*M[1][1]; 
int z = F[1][0]*M[0][0] + F[1][1]*M[1][0]; 
int w = F[1][0]*M[0][1] + F[1][1]*M[1][1]; 

F[0][0] = x; 
F[0][1] = y; 
F[1][0] = z; 
F[1][1] = w; 
} 

/*function that calculates F[][] raise to the power n and puts the 
result in F[][]*/
static void power(int F[][], int n) 
{ 
int i; 
int M[][] = new int[][]{{1,1},{1,0}}; 

// n - 1 times multiply the matrix to {{1,0},{0,1}} 
for (i = 2; i <= n; i++) 
	multiply(F, M); 
} 

This can be optimized to work in O(Logn) time complexity. We can do recursive multiplication to get power(M, n) in the previous method.

static int fib(int n) 
{ 
int F[][] = new int[][]{{1,1},{1,0}}; 
if (n == 0) 
    return 0; 
power(F, n-1); 
   
return F[0][0]; 
} 
   
static void multiply(int F[][], int M[][]) 
{ 
int x =  F[0][0]*M[0][0] + F[0][1]*M[1][0]; 
int y =  F[0][0]*M[0][1] + F[0][1]*M[1][1]; 
int z =  F[1][0]*M[0][0] + F[1][1]*M[1][0]; 
int w =  F[1][0]*M[0][1] + F[1][1]*M[1][1]; 
  
F[0][0] = x; 
F[0][1] = y; 
F[1][0] = z; 
F[1][1] = w; 
} 
   
static void power(int F[][], int n) 
{ 
if( n == 0 || n == 1) 
  return; 
int M[][] = new int[][]{{1,1},{1,0}}; 
   
power(F, n/2); 
multiply(F, F); 
   
if (n%2 != 0) 
   multiply(F, M); 
} 

Method 3 (O(log n) Time) Below is one more interesting recurrence formula that can be used to find nth Fibonacci Number in O(log n) time.

If n is even then k = n/2: F(n) = [2*F(k-1) + F(k)]*F(k)

If n is odd then k = (n + 1)/2 F(n) = F(k)*F(k) + F(k-1)*F(k-1) How does this formula work? The formula can be derived from the above matrix equation. fibonaccimatrix

Taking determinant on both sides, we get (-1)n = Fn+1Fn-1 – Fn2 Moreover, since AnAm = An+m for any square matrix A, the following identities can be derived (they are obtained from two different coefficients of the matrix product)

FmFn + Fm-1Fn-1 = Fm+n-1

By putting n = n+1,

FmFn+1 + Fm-1Fn = Fm+n

Putting m = n

F2n-1 = Fn2 + Fn-12

F2n = (Fn-1 + Fn+1)Fn = (2Fn-1 + Fn)Fn (Source: Wiki)

To get the formula to be proved, we simply need to do the following If n is even, we can put k = n/2 If n is odd, we can put k = (n+1)/2

public static int fib(int n) 
{ 
     
    if (n == 0) 
        return 0; 
          
    if (n == 1 || n == 2) 
        return (f[n] = 1); 
   
    // If fib(n) is already computed 
    if (f[n] != 0) 
        return f[n]; 
   
    int k = (n & 1) == 1? (n + 1) / 2 
                        : n / 2; 
   
    // Applyting above formula [See value 
    // n&1 is 1 if n is odd, else 0. 
    f[n] = (n & 1) == 1? (fib(k) * fib(k) +  
                    fib(k - 1) * fib(k - 1)) 
                   : (2 * fib(k - 1) + fib(k))  
                   * fib(k); 
   
    return f[n]; 
} 

Method 4 - Using a formula In this method, we directly implement the formula for the nth term in the Fibonacci series. Time O(1) Space O(1) Fn = {[(√5 + 1)/2] ^ n} / √5

static int fib(int n) { 
double phi = (1 + Math.sqrt(5)) / 2; 
return (int) Math.round(Math.pow(phi, n)  
                    / Math.sqrt(5)); 
} 

Reference: http://www.maths.surrey.ac.uk/hosted-sites/R.Knott/Fibonacci/fibFormula.html

We should first note that Fibonacci numbers (F(n)) grow very fast with n and cannot be represented in 64-bits for n larger than 93. So a program for computing them for such n needs to use additional mechanisms to operate on these large numbers. Now, considering only the count of (large-number) operations, the algorithm to sequentially compute them will require linear number of operations.

We can benefit from the below identity about Fibonacci numbers:

F(2m) = 2*F(m)*F(m+1) − (F(m))^2

F(2m+1) = (F(m))^2 + (F(m+1))^2

(a symbol like A^2 denotes square of A).

So, if we know F(m) and F(m+1), we can directly compute F(2m) and F(2m+1).

Consider the binary representation of n. Observe that starting with x = 1, we can make x = n by iteratively doubling and possibly adding 1 to x. This can be done by iterating over the bits of n, and checking if it is 0 or 1.

The idea is that, we can maintain F(x) in sync with x. In each such iteration, as we double x and possibly add 1 to x, we can also compute the new value of F(x) using the earlier value of F(x) and F(x+1), with above equations.

Since the number of iterations will be logarithmic in n, the total (large-number) operations are also logarithmic in n.