Nth Fibonacci number for n as big as 10^19?

Question

I am trying to make a program to find the nth Fibonacci number for 1 < n < 10^19.

Here is my code using dynamic programming.

memo = {}
def fib(n):
    if n in memo:
        return memo[n]
    if n <= 2:
        f = 1
    else:
        f = fib(n-1) + fib(n-2)
    memo[n]=f
    return f
print fib(input()) % 1000000007

My code does not seem to work for large numbers. I get invalid response error. Any suggestions?

Answer 1

Getting the Nth fibonacci number when N is 10^19 is not goign to work if you do it the naive way (at least i would guess it won't work).

There's a much better way to do it. And this technique works with lots of series' like this. It's called the Fibonacci Q Matrix.

Where

Think of it like this:

You have some matrix which transforms vector A into B:

Filling in those entries is easy. The special part is that this is now a matrix operator, and so if we want the 1000th Fibonacci number, we just need to do matrix multiplication.

You could do this with a loop, but it's going to take you quite a while to get all the way up to 10^19, and doing 10^19 matrix multiplications (even when they're small) is going to take a fair while too.

Instead, we take another shortcut. x^N can be rewritten as the product of power where they sum to N, i.e.

x**100 == x**90 * x**10

So the aim is to get large numbers in the indices without doing lots of calculations:

x**2 is just as difficult as x*x - they take the same amount of time. But x*x*x*x gives the same answer as (x**2)**2 while requiring an extra multiplication. The gains get more as you go to higher powers. So if you break down the exponent into powers of 2 (any power works, but this is the simplest case),

X**100 == X**64 * X**32 * X**4

i.e.

X**100 == (((((X**2)**2)**2)**2)**2)**2 + ...

So what you do, is work out the powers of two of the total power you want to reach, and then take the product of those powers of two of the Q matrix.

This seems to work for me:

fib_matrix = [[1,1],
              [1,0]]

def matrix_square(A, mod):
    return mat_mult(A,A,mod)


def mat_mult(A,B, mod):
  if mod is not None:
    return [[(A[0][0]*B[0][0] + A[0][1]*B[1][0])%mod, (A[0][0]*B[0][1] + A[0][1]*B[1][1])%mod],
            [(A[1][0]*B[0][0] + A[1][1]*B[1][0])%mod, (A[1][0]*B[0][1] + A[1][1]*B[1][1])%mod]]


def matrix_pow(M, power, mod):
    #Special definition for power=0:
    if power <= 0:
      return M

    powers =  list(reversed([True if i=="1" else False for i in bin(power)[2:]])) #Order is 1,2,4,8,16,...

    matrices = [None for _ in powers]
    matrices[0] = M

    for i in range(1,len(powers)):
        matrices[i] = matrix_square(matrices[i-1], mod)


    result = None

    for matrix, power in zip(matrices, powers):
        if power:
            if result is None:
                result = matrix
            else:
                result = mat_mult(result, matrix, mod)

    return result

print matrix_pow(fib_matrix, 10**19, 1000000007)[0][1]

And then, you can take it a step even further - it's just a 2x2 matrix, so we can diagonalise it, and then get the formula for the nth fibonacci number, just as a function of n - with no recursion. Like this:

As above, we compute the matrix that takes us from one step to the next:

And then the relationship to get from one set of numbers to the next:

where we can chain these matrix multiplications:

Where there's nothing to stop us going back all the way to the first fibonacci numbers:

now the game becomes "how do we raise that matrix to the power n" - which is exactly what's done in the code above. But there is a better way than the solution i pose above. We can decompose the Q-matrix into eigen values and vectors, a write it like so:

Where U is a unitary matricies that contain the eigen values of Q, and Λ is the matrix of corersponding eigenvalues. These eigenvalues and vecors are:

And then you use one of the standard advantages of this style of decomposition, where when you raise it to a power, the adjacent U matrix and it's inverse combine to give the unitary matrix, leaving you with a single U and it's inverse at the ends, with a chain of diagonal matrices in the middle, where raising these to a power is trivial:

So now we have all we need to write the nth Fibonacci number in terms of just a single formula, no recursion. I'll complete it tomorrow/some time later this week though...

Answer 2

At O(n) efficiency you'll never get there. Not specifically code-related, but Dijkstra's note "In honor of Fibonacci" describes a way to find F(n) in O(log(n)) efficiency.

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

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

That you could not only do, but still do recursively.

Answer 3

Python has a default recursion limit of 1000 (usually). To find out what the exact limit is on your system:

>>> import sys
>>> sys.getrecursionlimit()

Firstly, if you want to write this recursively and you're using Python 3.2 and above (which it doesn't look like you are, judging from the print statement) then you can use @functools.lru_cache(maxsize=128, typed=False) like so:

import functools

@functools.lru_cache()
def fib(n):
    if n <= 2:
        return 1
    else:
        return fib(n-1) + fib(n-2)

Having said that, this still won't be very fast for large numbers. The better way to do this is to write an iterative solution and all you need to "memoize", at any given time, is the last 2 numbers.

You can of course use the matrix form for even better performance.

Ultimately, for n being as large as 10**19 you're going to have a hard time writing anything that runs in Python without giving you an OverflowError.

Answer 4

I do not think you can go up to 1E19 with this, but here is how you avoid the double overflow and the recursion depth limit:

import decimal
import operator


def decimal_range(start, stop, step=1):
    """Provides an alternative to `xrange` for very high numbers."""
    proceed = operator.lt
    while proceed(start, stop):
        yield start
        start += step


def fib(n):
    """
    Computes Fibonacci numbers using decimal.Decimal for high 
    precision and without recursion
    """
    a, b = decimal.Decimal(0), decimal.Decimal(1)
    for i in decimal_range(0, n):
        a, b = b, a + b
    return a

On my machine, it took 26.5 s to compute 1E6, but I can't guarantee the correctness of the result:

In [26]: %time f2(n)
CPU times: user 26.4 s, sys: 130 ms, total: 26.5 s
Wall time: 26.5 s
Out[26]: Decimal('1.953282128707757731632014830E+208987')

The iterator is taken from this SO thread with minimal alterations, while the fib function can be found in this other thread.