nth fibonacci number in sublinear time

问题:

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

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

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

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

回答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

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

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

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

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

回答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

回答10:

see divide and conquer algorithm here

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

回答11:

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.

回答12:

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

回答13:

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.