nth fibonacci number in sublinear time

Question

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

Brose

One could argue that it's related to algorithms, since the OP makes a vague reference to algorithmic complexity... I'd still be curious what algorithm though. – Ardeb

The two answers below have the correct formula. On whether this question is programming-related: it's part of computer science. The apparatus used to derive the formula is known as "generating functions" and has an important role in algorithm analysis. – Casaleggio

@azheglov: While generating functions are useful, they are not needed to derive the closed-form expression for the Fibonacci sequence. – Plasmo

You have a problem that you want to solve for whatever reason, and you want to do it efficiently. Sometimes the required insight will be a new implementation, sometime an algorithm, and sometimes mathematics. There is no need to decry the situation as "not programming related" every time the latter happens. – Byrd

The size of the result is linear in n. Therefore there is no such algorithm. Of course that doesn't invalidate any of the nice answers below that compute Fibonacci numbers using O(log n) arithmetic operations. – Swats

@Swats Yes. Because they are computing approximations. – Phthisis

It looks like #30596344 – Minnick

It depends what "time" is: if "time" is measured in arithmetic operations, then you can compute the result in log(n) time. If "time" is measured in bit cost, then you can't. Note that in much of complexity theory, the choice of "time" is flexible (for example: sorting algorithms use "time" to mean pairwise comparisons, and also hash table analysis use "time" to mean comparisons). Neither of these accurately describe time as performed on a "real" computer (because computing a hash is not O(1) for arbitrary large objects), but no one complains about them). – Ewer

Answer 1

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);
}

Plasmo

As an aside, I'd probably advocate making those variables const rather than static, just to be safe. – Ardeb

@Matthew Scharley: Yes, they probably should, but then I would have to key in the decimal value for 1 / Math.Sqrt(5) and (1 + Math.Sqrt(5)) / 2 as C# won't let const values be computed using Math.Sqrt. For our purposes here, the formula are clearer for understanding than the mysterious decimal values would be. – Plasmo

@Jason: that's what comments are for. But as it's static the difference is so small it's not worth messing with. – Walkup

@Joel Coehoorn: I agree with you. I was trying to keep it simple here. – Plasmo

@Json I have not downvoted you, but others may be doing so because your answer suggests that the Nth fibonacci number can be computed in O(log n) time, which is false. Your code is computing an approximation. Your code would be at least O(n) in arbitrary precision, because the length of the answer is O(n). – Phthisis

@PeterAllenWebb: The formula provided is not an approximation. The nth Fibonacci number is equal to the floor of phi^n / sqrt(5) + 1/2 where phi = (1 + sqrt(5)) / 2. This is a fact. Second, I understand the point that others are making about the length of the answer being O(n) but I have added a remark to my answer assuming that the primitive mathematical operations take constant time (I know they are not unless you bound the inputs). My point is that we can find the nth Fibonacci number in O(log n) arithmetic operations. – Plasmo

@Jason: Assuming that exponentiation is O(1) too makes the whole algorithm O(1). That would be nice, however, exponentiation is not O(1) and neither are the other primitive mathematical operations. So in short, the formula is nice, but it doesn't calculate the result in sub-linear time. – Dyak

Here's something interesting: double foo = 1.6; printf ("%.40f\n", foo); Prints: 1.6000000000000000888178419700125232338905 Now this is to be expected because of IEEE 754 floating point. Let's try your constants above: double foo = 1.6180339887498948482045868343656; printf ("%.40f\n", foo); Prints: 1.6180339887498949025257388711906969547272 So, the question is when does this rounding error start to matter? – Ancestor

@Jason: The formula is not an approximation, but the code is an approximation (except on an imaginary C# implementation in which Math.Pow(…) has infinite precision, in which case the code is O(n)). – Byrd

@ShreevatsaR: Actually, note the use of Math.Floor. The result of executing this code, even on a pathetic machine that doesn't have infinite precision, is the nth Fibonacci number. – Plasmo

@Jason: Nope. Run your code on n=1000 (for which the Fibonacci number 43466...849228875 has a measly 209 digits) and tell me if you get all digits right. For Math.Floor to get the integer part right, those many digits have to be accurately calculated by Math.Pow. In fact, on my C++ implementation, even the 16-digit F_{74} = 130496954492865 is calculated incorrectly, even though the integer 130496954492865 can be represented exactly (with long long), and I'd be surprised if C# gets much more digits than that. – Byrd

@Jason's code is correct for the first 46 Fibonacci numbers then it overflows. The conventional a,b = b, a+b algorithm does the same. However if you change int to long in the C#, then it makes a mistake at n=71. – Tc

@Plasmo it is giving wrong results for higher values – Stanislaw

@Byrd Jason's precise implementation happens to be approximate, but it is relatively easy to implement this formula without floating point arithmetic at all. See for example this answer – Batholomew

@Scott Sure, if you have arbitrary-precision integers, as in your Python answer. That was exactly my point (see comment of "Oct 24 '11 at 3:07") -- to repeat, the formula is exact, but the code in this answer (the C# implementation) is not. – Byrd

(for most computer language implementations) you can make the code slightly more accurate by replacing 0.5 with 0.28 | 0.27 doesn't make it more accurate | 0.29 doesn't make it more accurate – Larimore

Answer 2

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   Φ ]

       = ( √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 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.

Kenyatta

+1 - Awesome stuff, as usual. What did you use to typeset it? LaTeX? – Affusion

It is copy-pasted from the book of Algebra of Gilbert Strang, or from other good book of Linear Algebra. – Monoclinous

@Monoclinous it was not 'copy pasted', but done as an exercise to check I could still remember my lin a, with some reference to Open University course notes and wikipedia. – Kenyatta

I took the course of L Algebra with Gilbert Strang, and there it was identical. Quite so, the problem of expressing recursion via matrix decomposition is classical, and can be found in any good text book / course. – Monoclinous

This is a very long time later but saying "this gives an O(1) method of raising M to the nth power" is wrong. The Eigenvalue decomposition method may reduce the number of multiplications required to get a result but raising a number to a power n remains a O(log(n)) operation (when assuming * is O(1)). – Khichabia

@Khichabia yes, I mention that an exact method is O(log(n)) but if you want an approximate method you can use eigenvalue decomposition and then floating point math is O(1) – Kenyatta

Oh I see, simple floating-point exponantiation to get the most significant few digits, is it? Just wondering how far this method can go before the first out-of-sequence number? – Khichabia

Answer 3

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.

Dyak

Can you explain this further? 2^n has n+1 significant digits but can be generated in constant time. – Sanative

@Sumit: That depends on what output you want. If, for example, you want to store all of your number's digits, then even for 2^n you will have to allocate and clear a buffer of the appropriate size and that would take linear time. – Dyak

@yairchu: Let me rephrase this, if I understand it correctly. In theory calculating fib_n requires calculating n digits so for any arbitrary n it will take O(n) time. However, if fib_n < sizeof(long long) then we can calculate fib_n in O(log n) time since the machine architecture is providing a parallel mechanism of setting the bits. (For example, int i = -1; requires setting 32-bits but on a 32-bit machine all the 32 bits can be set in constant time. – Sanative

@Sumit: If you only want to support results that fit in 32-bits, then you can also have a lookup table for these first 48 results of the series. That's obviously O(1), but: Doing big-O analysis for a bounded N is silly, as you can always incorporate anything into the constant factor. So my answer refers to unbounded input. – Dyak

@yairchu: Could you demonstrate your logic for a well-known example such as O(n*log n) for comparison-based sorting of a sequence of n numbers where each number has O(log n) digits? – Gunlock

This is right or wrong depending on what you intend "time" to mean. For sorting (or hash table look ups), "time" means the number of comparisons. In the question it could mean arithmetical operations. In this answer it's taken to mean something like digit-wise operations. – Ewer

That depends on your choice of basis. Computing the Fibonacci numbers in base Phi(or even more trivially, in a fibonacci encoding) should be O(1). – Theriault

@saolof: Do you often use non-integer bases? – Dyak

Often may not be the best way to describe it, but base phi is the third most common basis I use, after Decimal and binary (counting octal and hex as binary). In some cases, you want to use a basis as close to 1 as possible. But the smallest integer larger than 1 is 2. If that's too big, you have to use a non-integer basis and for these cases the golden ratio is the most common choice since integers have finite representations. – Theriault

@Theriault when would you like to use a basis as close to 1 as possible? and wouldn't integers have finite representation in base sqrt(2), which is even closer to 1, as well? how is phi any better than sqrt(2)? – Dyak

Integers will indeed have a finite representation in base sqrt(2), but it will just be zero on odd digits, i.e. equivalent to base 2. If any of the odd digits in base sqrt(2) is nonzero, you have an irrational number. One case where you may want base phi is in ADCs when converting continuous signals to analog. Afaik this is the "industrial" application of base phi, where it is used to reduce coarse graining when rounding the signal. Personally though, I used base phi and fibonacci encodings as a notationally convenient way to work with Fibonacci anyon representations of the braid group. – Theriault

@Theriault Very interesting (en.wikipedia.org/wiki/Beta_encoder). Still haven't quite learned it but it's cool to know that such things have had practical use. – Dyak

@jfs: it's meaningful to talk about complexity classes for sorting an unbounded number of fixed-size objects. Or even to talk about a sorting a bounded number of fixed-size objects, because it grows slowly enough for asymptotic effects to matter on practical computers. – Apodaca

@jfs: But Fib(n) grows exponentially and thus hits any fixed width so quickly that we need to consider implementation details like CPU microarchitecture and constant factors, not just complexity class. e.g. Indexed branch overhead on X86 64 bit mode, and Assembly Language (x86): How to create a loop to calculate Fibonacci sequence. Also related: Write the fastest Fibonacci Fib(20mil) using arbitrary-precision techniques (with O(log2(n)) extended-precision steps using GMP). – Apodaca

Answer 4

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
             p ← p² + q²
             q ← 2pq + q²
             count ← count ÷ 2
        Else
             a ← bq + aq + ap
             b ← bp + aq
             count ← count - 1
        End If
    End While

    Return b
End Function

Rudy

here's an implementation in Python (to be used with twisted framework). – Gunlock

"If Even(count) Then" should be "If Odd(count) Then" – Zymo

@MonirulIslamMilon if even(count) is correct. The sequence starts with zero (zeroth Fibonacci number is zero): 0,1,1,2,3,5,8,13,... – Gunlock

The book link is now: mitpress.mit.edu/sites/default/files/sicp/full-text/book/… – Pokeberry

Late comment, but the variables p and a are overwritten before being used to calculate q and b. To avoid this issue, pre-calculate terms and change order of p and q assignments : | qq = q·q | q = 2·p·q + qq | p = p·p + qq | ... | aq = a·q | a = b·q + aq + a·p | b = b·p + aq | . – Stateroom

Answer 5

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.

Liselisetta

Use the eigen decomposition of M to calculate M^n efficiently. – Kenyatta

Proposed method is fine for calculations in integers (probably with long arithmetic). Approach with eigen decomposition is not interesting: if you don't need integer calculations, then use formula from Jason's answer. – Headrail

@Konstantin The formula from Jason's answer is the result given by eigen decomposition, so you're contradicting yourself. – Kenyatta

@Pete Kirkham That formula can be obtained by several methods: characteristics equation, eigen decomposition, proof by induction. I'm not sure, that eigen decomposition is the easiest one. In any case it is well-known, and it is easier to use it immediately – Headrail

Answer 6

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;
    }

Chisel

You can avoid the need to compute to two exponentials by using the fact that |1 - phi|^n / sqrt(5) < 1/2 when n is a nonnegative integer. – Plasmo

Didn't know that adjustment always have used the other form, but that is a nice optimisation – Chisel

Approximation of the result the correct solution involves matrix multiplication. – Atonsah

Answer 7

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.

Humeral

This code isn't so efficient. Imagine that the fibs[] array is only size 10 and you call Fib(101). Fib(101) calls Fib(51) and Fib(50). Fib(51) calls Fib(26) and Fib(25). Fib(50) calls Fib(25) and Fib(24). So Fib(25) was called twice, which is a waste. Even with fibs up to 500, you will have the same problem with Fib(100000). – Ida

Answer 8

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.

Ida

I want to stress here that you want to compute F(2n) and F(2n+1) in function of F(n) and F(n-1). You did not indicate what you want to do. – Monoclinous

Answer 9

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

Sura

Answer 10

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.

Tc

Answer 11

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+p1*phi)(q0+q1*phi) = p0q0 + (p0q1+q1p0)phi + p1q1*phi^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.

Ewer

Answer 12

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)

Magnesite

This is in linear time though, the question specifically asks how to achieve sublinear time (which is possible using a sort of closed-form solution). – Numbskull

Answer 13

see divide and conquer algorithm here

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

Totalizer

bad link - returns 403 forbidden response – Fumikofumitory

Answer 14

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

Batholomew

Answer 15

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

Seaworthy

Answer 16

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.

Chimp

How many of the pre-existing answers even to this question mentioned this same method? The question asked for sub-linear time and you argued about large-number operations - What is the asymptotic time complexity for a RAM? See Accipitridae's comment, too. – Trailer

Answer 17

c++14 sample implementation:

uint64_t fib(unsigned n) {
    static struct FIB_NUMBERS {
        uint64_t numbers[100] {0, 1};
        unsigned max_n = 2;
        constexpr FIB_NUMBERS() {
            for (; max_n<100; ++max_n) {
                numbers[max_n] = numbers[max_n-1] + numbers[max_n-2];
                if (numbers[max_n] < numbers[max_n-1]) {
                    // overflow
                    numbers[max_n--] = 0;
                    break;
                }
            }
        }
    } fib;

    if (n <= fib.max_n)
        return fib.numbers[n];
    return 0;
}

Compiler will generate all fib numbers until n = 93 Runtime is just a lookup. To detect overflow from the caller:

auto res = fib(n);
if (res==0 && n>0) overflow_msg();

Buster