Computing huge Fibonacci numbers

In this post, we will show one rather efficient technique for computing large Fibonacci numbers. The resulting Java implementation computes F_{1 000 000} in less than two seconds on a 2.5 GHz CPU. The key idea is that

\displaystyle  A^n =   \begin{pmatrix}  1 & 1 \\  1 & 0  \end{pmatrix}^n =  \begin{pmatrix}  F_{n + 1} & F_n \\  F_n       & F_{n - 1}  \end{pmatrix},

where

\displaystyle  F_n =   \begin{cases}  0 & \text{if } n = 0, \\  1 & \text{if } n = 1, \\  F_{n - 1} + F_{n - 2} & \text{if } n > 1.  \end{cases}

Let us use induction in order to prove the matrix equation.

The base case

\displaystyle  \begin{pmatrix}  1 & 1 \\  1 & 0  \end{pmatrix}^1  =  \begin{pmatrix}  1 & 1 \\  1 & 0  \end{pmatrix}  =   \begin{pmatrix}   F_2 & F_1 \\  F_1 & F_0  \end{pmatrix}.

Inductive step

Suppose that

\displaystyle  \begin{pmatrix}  1 & 1 \\  1 & 0  \end{pmatrix}^n  =  \begin{pmatrix}  F_{n + 1} & F_n \\  F_n       & F_{n - 1}  \end{pmatrix}.

Now we have that

\displaystyle  \begin{aligned}  \begin{pmatrix}  1 & 1 \\  1 & 0   \end{pmatrix}^{n + 1}   &=  \begin{pmatrix}  1 & 1 \\  1 & 0  \end{pmatrix}^n  \begin{pmatrix}  1 & 1 \\  1 & 0  \end{pmatrix} \\  &=  \begin{pmatrix}  F_{n + 1} & F_n \\  F_n       & F_{n - 1}  \end{pmatrix}  \begin{pmatrix}  1 & 1 \\  1 & 0  \end{pmatrix} \\  &=  \begin{pmatrix}  F_{n + 1} + F_n & F_{n + 1} \\  F_n + F_{n - 1} & F_n  \end{pmatrix} \\  &=  \begin{pmatrix}  F_{n + 2} & F_{n + 1} \\  F_{n + 1} & F_n  \end{pmatrix},  \end{aligned}

which concludes the proof. Since we want to deal with large Fibonacci numbers, we have to opt to java.math.BigInteger. Using BigInteger not just allows dealing with integers that do not fit in primitive integer types, but also provides more efficient methods for multiplying two integers: as of now, JDK 8 implementation of BigInteger relies on two distinct multiplication algorithms

(The naïve multiplication runs in \Theta(n^2).) Another optimization is to compute the matrix exponent by divide-and-conquer technique that requires \Theta(\log n) actual multiplications: basically, if n is even, compute recursively B = A^{n / 2} and return B^2; if n is odd, compute recursively C = A^{n - 1} and return AC. (The recursion is terminated at n = 1, where we return the input matrix unmodified.) Since F_n = \Theta(c^n) for some c > 1, the length of the integer is \Theta(\log F_n) = \Theta(\log c^n) = \Theta(n). Also, the worst case complexity of multiplying two BigIntegers of length n is \mathcal{O}(n^{1.585}). As we need logarithmic amount of matrix multiplications and each matrix multiplication involves a constant number of BigInteger multiplications, we have the upper bound of \mathcal{O}(n^{1.585} \log n) for computing F_n.

Source code

import java.math.BigInteger;

public class LargeFibonacciNumbers {

    public static BigInteger fibonacci(int n) {
        if (n < 0) {
            throw new IllegalArgumentException(
                    "The (" + n + ")th Fibonacci number is not defined.");
        }

        switch (n) {
            case 0:
                return BigInteger.ZERO;

            case 1:
                return BigInteger.ONE;
        }

        BigInteger[][] matrix = new BigInteger[2][2];

        matrix[0][0] = BigInteger.ONE;
        matrix[0][1] = BigInteger.ONE;
        matrix[1][0] = BigInteger.ONE;
        matrix[1][1] = BigInteger.ZERO;

        return fibonacciMatrix(matrix, n)[0][1];
    }

    private static BigInteger[][] multiply(BigInteger[][] matrix1, 
                                           BigInteger[][] matrix2) {
        BigInteger[][] ret = new BigInteger[2][2];

        ret[0][0] = matrix1[0][0].multiply(matrix2[0][0])
               .add(matrix1[0][1].multiply(matrix2[1][0]));

        ret[0][1] = matrix1[0][0].multiply(matrix2[0][1])
               .add(matrix1[0][1].multiply(matrix2[1][1]));

        ret[1][0] = matrix1[1][0].multiply(matrix2[0][0])
               .add(matrix1[1][1].multiply(matrix2[1][0]));

        ret[1][1] = matrix1[1][0].multiply(matrix2[0][1])
               .add(matrix1[1][1].multiply(matrix2[1][1]));

        return ret;
    }

    private static BigInteger[][] 
        fibonacciMatrix(BigInteger[][] matrix, int n) {
        if (n == 1) {
            // End the recursion.
            return matrix;
        }

        if ((n & 1) == 1) {
            // 'n' is odd.
            return multiply(matrix, fibonacciMatrix(matrix, n - 1));
        } 

        // 'n' is even.
        BigInteger[][] tmp = fibonacciMatrix(matrix, n >> 1);
        return multiply(tmp, tmp);
    }

    public static void main(String[] args) {
        int n = -1;

        if (args.length > 0) {
            try {
                n = Integer.parseInt(args[0]);
            } catch (NumberFormatException ex) {
                System.err.println("Could not parse \"" + args[0] + 
                                   "\" as an integer.");
                return;
            }
            
            try {
                System.out.println(fibonacci(n));
            } catch (IllegalArgumentException ex) {
                System.err.println(ex.getMessage());
            }
        } else {
            System.out.println("Usage: java -jar File.jar N");
            System.out.println("Where N is the index of the Fibonacci number " +
                               "to compute.");
        }
    }
}
Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s