The Nth Fibonacci Number in O(log N)

Comments

1. 

   jannschu 3:45AM on Wednesday, Dec 17, 2014

   Computing the n-th power of the diagonalized matrix would involve non-integer arithmetic, thus rounding errors. One could argue that the error should not affect the integer part, but I’d say that depends on the size of n.

   The power function is also not O(1) in n. You would get the same runtime O(log(n)) with a clever implementation.

   Reply
2. 

   Árni Steingrímur Sigurðsson 4:57AM on Wednesday, Dec 17, 2014

   Bottom of def get_number renders as the Registered Trademark unicode instead of ®.

   You can be a little more concise with the code by modifying how __memo is set up.

   def __init__(self):
           self.__memo = {0:0,1:1}
   

   Then you can do:

   
   def __get_matrix_power(self, M, p):
           """Matrix exponentiation (it is expected that p that is equal to the power of 2)."""
           return self.__memo.get(p,self.non_memoized(M,p))
   
   def non_memoized(self, M,p):
           K = self.__get_matrix_power(M, int(p/2))
           R = self.__multiply_matrices(K, K)
           self.__memo[p] = R
           return R
   
   

   I know it essentially does the same thing, IMHO it’s a little cleaner code thou ;)

   Reply
3. 

   Árni Steingrímur Sigurðsson 5:09AM on Wednesday, Dec 17, 2014

   Also, you should memoize the matrix multiplication, that’s a heavy step, and one for repeated invocations would save a ton of time.

   Reply
4. 

   Randall Mason 2:45AM on Thursday, Dec 18, 2014

   Why do you memoize the matrix method and not the iterative method? When I memoize the iterative method, I end up with about constant time for each of the steps: Around 1000 gives Iterative: 0.054137468338 Matrix: 1.92105817795

   This all stays in memory for my machine.

   You can find the improved fork of your code here: gist.github.com/ClashTheBunny/762a497d920a24aada69

   Reply
5. 

   distorted-mdw 10:03AM on Tuesday, Jan 6, 2015

   Because the exponentiation method lets you compute an arbitrary element of the Fibonacci sequence without calculating most of the ones before it. Your `benchmark’ calculates them in ascending order, so obviously a straightforward memoizing/iterative implementation does well. Instead, it should just calculate the millionth Fibonacci number (say), ab initio.

   Reply
6. 

   Juan Lopes 9:41AM on Friday, Dec 19, 2014

   Thinking about it, I realized that even the implementation described here isn’t really O(log n). Here’s why: it’s rather easy to prove that the nth term of the Fibonacci sequence has O(n) digits. So, even if you are performing only O(log n) multiplications, each multiplication is actually O(n), so the algorithm is actually O(n log n).

   I made a benchmark that shows that in practice:

   docs.google.com/spreadsheets/d/1mDIuAUrI6lzZxDf2HEqHbaYjddfK0hijd_r9qBOOEoQ

   Reply
7. 

   Kukuruku Hub 2:13PM on Friday, Dec 19, 2014

   You should write an article — «Response to the Nth Fibonacci Number in O(log N)»! IT would be interesting to read.

   Reply
8. 

   distorted-mdw 10:00AM on Tuesday, Jan 6, 2015

   You’re doing too much work in the matrix multiplication: the matrices are always symmetrical, but you’re not taking advantage of that. But there’s a slicker way, which involves some slightly fancier mathematics.

   Suppose that we know that t^2 = t + 1 for some t. (Equivalently, we’ll work in the quotient ring R = Z[t]/(t^2 — t — 1).) Now, let’s look at the powers of t. By using the rule that t^2 = t + 1, we will end up with only units and coefficients of t. (That this is the characteristic polynomial of the matrix discussed in the main article is, of course, not a coincidence.)

   t^0 = 1 t^1 = t t^2 = t + 1 t^3 = t^2 + t = (t + 1) + t = 2 t + 1 t^4 = 2 t^2 + t = 2 (t + 1) + t = 3 t + 2 t^5 = 3 t^2 + 2 t = 3 (t + 1) + 2 t = 5 t + 3

   Hmm. There’s a pattern: t (a t + b) = a t^2 + b t = a (t + 1) + b t = (a + b) t + a. I feel a theorem coming on: t^k = F_k t + F_{k-1} for all integers k. Proof. (a) t^1 = t = F_1 t + F_0. (b) Suppose the claim holds for k; then t^{k+1} = t t^k = (F_k + F_{k-1}) t + F_k = F_{k+1} t + F_k. © Notice that t (t — 1) = t^2 — t = (t + 1) — t = 1, so t^{-1} = t — 1; then, supposing again that the claim holds for k, we have t^{k-1} = (t — 1) (F_k t + F_{k-1}) = F^k t^2 + F_{k-1} t — F_k t — F_{k-1} = F_k (t + 1) + F_{k-1} t — F_k t — F_{k-1} = F_{k-1} t + (F_k — F_{k-1}) = F_{k-1} t + F_{k-2}, as required. []

   Now, we can compute powers of t using the ordinary square-and-multiply method (or fancier methods involving windowing, combs, or whatever). The important pieces are:

   (a t + b)^2 = a^2 (t + 1) + 2 a b t + b^2 = (2 a b + a^2) t + (a^2 + b^2) t (a t + b) (c t + d) = a c (t + 1) + (a d + b c) t + b d = (a d + a c + b c) t + (a c + b d)

   Here’s a simple left-to-right exponentiation algorithm. Suppose we want to compute x^n, and we know that n = 2^k u + v, for 0 <= v < 2^k, and x^u = y. (Initially, we can arrange this with 2^k > n, so u = 0, v = n, and x = 1.) Now, suppose that v = 2^{k-1} b + v’, where b is either zero or one, and v’ < 2^{k-1}. Let u’ = 2 u + b. Then n = 2^{k-1} u’ + v’, and x^{u’} = x^{2u+b} = y^2 x^b. Continue until k = 0.

   Annoyingly, this means we need to count the number of bits in n, i.e., k = L(n), with 2^{i-1} <= n < 2^i. We can start with an upper bound, finding the smallest 2^{2^j} > n. Suppose, inductively, that we know how to determine L(m) for 0 <= m < 2^i, and that n < 2^{2i}. If n < 2^i, then we’re done; otherwise write n = 2^i n’ + r where 0 <= r < 2^i. Since n < 2^{2i}, we know that n’ <= n/2^i < 2^i, so we can determine k’ = L(n’). Let k = k’ + i: I claim that 2^{k-1} <= n = 2^i n’ + r < 2^k. Taking the lower bound first, we’re given that 2^{k’-1} <= n’, so 2^{k-1} = 2^i 2^{k’-1} <= 2^i n’ <= 2^i n’ + r. But also, n’ < 2^{k’}, so n’ + 1 <= 2^{k’} and 2^i n’ + 2^i <= 2^i 2^{k’} = 2^k. But r < 2^i, so n = 2^i n’ + r < 2^k as required. []

   Enough theory! Now let’s turn this into actual code!

   
   #! /usr/bin/python
   
   def power(x, n, id, sqr, mul):
     """
     Compute and return X^N for nonnegative N.
   
     This is an abstract operation, assuming only that X is an element of a
     monoid: ID is the multiplicative identity, and the necessary operations are
     provided as functions: SQR(Y) returns Y^2, and MUL(Y, Z) returns the
     product Y Z.
     """
   
     ## First, we need to find a 2^{2^i} > n.
     i = 1
     while (n >> i) > 0: i = 2*i
   
     ## Refine our estimate, so that 2^{k-1} <= n < 2^k.
     k, m = 1, n
     while m > 1:
       i = i >> 1
       mm = m >> i
       if mm > 0: m, k = mm, k + i
   
     ## Now do the square-and-multiply thing.
     y = id
     for i in xrange(k - 1, -1, -1):
       y = sqr(y)
       if (n >> i)%2: y = mul(x, y)
   
     ## We're done.
     return y
   
   def fibonacci(n):
     """
     Compute and return the Nth Fibonacci number.
   
     Let F_0 = 0, F_1 = 1, F_{k+1} = F_k + F_{k-1}, for all integers k.  Note
     that this is well-defined even for negative k.  We return F_N.
     """
   
     ## Auxiliary functions for working in our polynomial ring.
     def poly_sqr((a, b)):
       a2 = a*a
       return 2*a*b + a2, a2 + b*b
     def poly_mul((a, b), (c, d)):
       ac = a*c
       return a*d + b*c + ac, ac + b*d
   
     ## Do the job.  For negative indices, we take powers of t^{-1}.
     if n < 0: return power((1, -1), -n, (0, 1), poly_sqr, poly_mul)
     else: return power((1, 0), n, (0, 1), poly_sqr, poly_mul)
   
   ## Command-line interface.
   if __name__ == '__main__':
     import sys
     try: n = (lambda _, n: int(n))(*sys.argv)
     except TypeError, ValueError:
       print >>sys.stderr, 'Usage: fib N'
       sys.exit(1)
     print fibonacci(n)[0]
   
   

   Reply
9. 

   Danny Hermes 3:09AM on Wednesday, Dec 17, 2014

   Just diagonalize the matrix and do it in constant time.

   Reply
10. 

    George van den Driessche 1:28PM on Wednesday, Dec 17, 2014

    If you do the matrix power more cunningly, you don’t need to memoize anything. Just use this:

    
      def mat_pow(m, e):
        if e == 1:
          return m
        else:
          m2 = mat_pow(m, e/2)
          if (e % 1) == 0:
            return m2 * m2
           else:
            return m2 * m2 * m
    
    

    Reply
11. 

    Juan Lopes 8:44AM on Friday, Dec 19, 2014

    
    def mul(A, B):
        return (A[0]*B[0] + A[1]*B[2], A[0]*B[1] + A[1]*B[3], 
                A[2]*B[0] + A[3]*B[2], A[2]*B[1] + A[3]*B[3])
        
    def fib(n):
        if n==1: return (1, 1, 1, 0)
        A = fib(n>>1)
        A = mul(A, A)
        if n&1: A = mul(A, fib(1))
        return A
    
    

    *cap obvious flies away*

    Reply
12. 

    reneleyva 2:40PM on Thursday, Apr 23, 2015

    I really enjoyed this article, thanks!

    Reply
13. 

    Coding Skull 7:30PM on Friday, Dec 11, 2015

    You can also use the fact that, if f(n) is the nth Fibonacci number (using f(0) = 1), then f(n + k) = f(n — 1)f(k — 1) + f(n)f(k)

    This is not hard to see by looking at the sequence, assuming that a = f(n), b = f(n + 1) for some n: a, b, a + b, a + 2b, 2a + 3b, 3a + 5b, 5a + 8b, 8a + 13b, 13a + 21b,…

    You can use this to do a similar ‘exponentiation by squaring’ thing:

    def add((a, b, c), (x, y, z)):
        v = a * y + b * z
        w = b * y + c * z
        return (w - v, v, w)
    
    def fib2(n):
        (r, p) = ((1, 1, 2), (0, 1, 1))
        while not(n == 0):
            if n % 2: r = add(r, p)
            p = add(p, p)
            n = n >> 1
        return r[0]
    

    Which is logarithmic in n if you assume that multiplication, addition, etc. are constant time (which is only true for word-sized integers, so in general this is false).

    Reply
14. 

    Angelo Catalani 6:56AM on Thursday, Sep 8, 2016

    However in my opinion time complexity is O(z) since it’s still linear, In fact ,we consider n as a integer so that its real input size is 2^z .Finally O(log(n)) can be written as O(z);

    Reply
15. 

    Angelo Catalani 7:03AM on Thursday, Sep 8, 2016

    The matrix version is better than the iterative version because the latter is 0(n) and so O(2^z)

    Reply