nth fibonacci number in sublinear time

Question

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

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

Solution 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 

Λ

n

Λ

n

M

n

Λ

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.

Solution 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.

Solution 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

Solution 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.

View more solutions

Comments

[Matthew Scharley]

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.

[azheglov]

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.

[jason]

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

[sidgeon smythe]

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.

[Accipitridae]

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.

[PeterAllenWebb]

@Accipitridae Yes. Because they are computing approximations.

[Rodolfo]

It looks like stackoverflow.com/questions/30595844/…

[Paul Hankin]

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).

[jason]

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.

[Matthew Scharley]

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

[JDunkerley]

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

[jason]

@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.

[Pete Kirkham]

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

[Joel Coehoorn]

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

[jason]

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

[Konstantin Tenzin]

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.

[Pete Kirkham]

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

[Konstantin Tenzin]

@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

[duffymo]

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

[PeterAllenWebb]

@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).

[jason]

@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.

[yairchu]

@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.

[Hank]

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?

[sidgeon smythe]

@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)).

[jason]

@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.

[sidgeon smythe]

@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.

[Sumit]

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

[yairchu]

@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.

[Sumit]

@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.

[yairchu]

@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.

[Colonel Panic]

@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.

[jfs]

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

[Eyal]

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).

[cerkiewny]

Approximation of the result the correct solution involves matrix multiplication.

[alinsoar]

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

[Pete Kirkham]

@alinsoar 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.

[alinsoar]

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.

[alinsoar]

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.

[jfs]

@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?

[Md. Monirul Islam]

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

[Paul Hankin]

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.

[jfs]

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

[saolof]

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).

[yairchu]

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

[saolof]

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.

[yairchu]

@saolof 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)?

[saolof]

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.

[yairchu]

@saolof 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.

[deenbandhu]

@jason it is giving wrong results for higher values

[Peter Cordes]

@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.

[Peter Cordes]

@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).

[Lee D]

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

[Romeo Valentin]

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).

[Him]

@ShreevatsaR 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

[sidgeon smythe]

@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.

[rcgldr]

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 | .

[beauxq]

(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

[greybeard]

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.