Fibonacci Matrix

I will now discuss an efficient implementation for generating the fibonacci sequence.

If you are not familiar with matrix multiplication, you can read the next two sections (algorithm and analysis) and skip the rest. You can come back to read the code after reading the Linear Algebra chapter.

\[ \begin{pmatrix} F_{n+1} \\ F_{n} \end{pmatrix} = \begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix} \begin{pmatrix} F_{n} \\ F_{n-1} \end{pmatrix} \\ \]

The above equation can be verified trivially. And by applying the equation to itself repeatedly, we end up with:

\[ \begin{equation} \begin{pmatrix} F_{n+1} \\ F_{n} \end{pmatrix} = \begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix}^n \begin{pmatrix} 1 \\ 0 \end{pmatrix} \end{equation} \]

Exponentiation algorithm

Multiplications are expensive. Our ability to compute fibonacci numbers quickly is predicated on computing powers fast. I.e finding the value of \(A^n\) where

\[ A = \begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix} \]

This algorithm is quite intuitive. Suppose you want to calculate \(3^{15}\). You either have to multiply 3 by itself 15 times. Or do the following:

\[ \begin{align*} 3^{15} &= 3 \times 3^{14} \\ &= 3 \times (3^7)^2 \\ &= 3 \times \left(3 \times (3^3)^2\right)^2 \\ &= 3 \times \left(3 \times (3 \times 3^2)^2\right)^2 \end{align*} \]

We can start squaring and reduce this value from inside.

\[ \begin{align*} 3^{15} &= 3 \times \left(3 \times (3 \times 9)^2\right)^2 \\ &= 3 \times \left(3 \times 27^2\right)^2 \\ &= 3 \times \left(2187\right)^2 \\ &= 3 \times 4782969 \\ &= 14348907 \end{align*} \]

We have done only 5 multiplications instead of 15 by repeated squaring! Let's generalize this algorithm.

\[ \begin{equation*} A^n = \begin{cases} \left(A^2\right)^{\frac{n}{2}} & \text{if}\ \text{n is even} \\ A \times \left(A^2\right)^\frac{n-1}{2} & \text{if}\ \text{n is odd} \end{cases} \end{equation*} \]

Big picture: we compute \(A^2\) and pass it recursively to the pow function with n halved.

Analysis

From the above observations, we can write the time complexity equation:

\[ \begin{align*} T(n) &= C + T\left(\frac{n}{2}\right) \\ &= C + C + T\left(\frac{n}{4}\right) \\ &= kC + T\left(\frac{n}{2^k}\right) \end{align*}\]

When \(n = 2^k\) or \(\log_{2}(n) = k\) we reach the base case.

\[ \begin{align*} T(n) &= C\log_2(n) + T(1) \\ T(n) &= \Theta(\log(n)) \end{align*}\]

This is a very fast logarithmic algorithm that can potentially compute the trillionth fibonacci number in under a second. The memory usage is once again constant. We declare some values and keep updating them.

Finally getting into code...

Equation (1) when translated into Odin looks like this:

fibonacci_matrix :: proc(n: int) -> int {
    f := matrix[2, 1]int{
        0,
        1,
    }

    a := matrix[2, 2]int{
        1, 1,
        1, 0,
    }

    result := matrix_pow(a, n) * f
    return result[0, 0]
}

Most of it should be self explanatory. You can also specify the type on the left of = during assignment.

f : matrix[2, 1]int = {
    0, 
    1
}

result[0, 0] returns the first element of the column matrix. We just need to implement matrix_pow that uses the above algorithm.

matrix_pow :: proc(m: matrix[2, 2]int, pow: int) -> matrix[2, 2]int {
    result := 1
    n := pow
    m := m

    for n > 1 {
        if n % 2 != 0 {
            result = m * result
            n -= 1
        }
        m = m * m
        n /= 2
    }

    return result * m
}

We could have initialized result to the identity matrix. But it's not necessary. Odin supports multiplying a matrices with a number (scalars).

result := matrix[2, 2]int{
        1, 0,
        0, 1,
}

for n > 1 says that the loop has to continue as long as the condition is satisfied. Odin doesn't have while. % is the remainder operator. n % 2 == 0 checks if n is odd.

n -= 1 is an shorter syntax for n = n - 1. We are subtracting 1 from n and updating the value of n. Remember = is mutation not mathematical equality. Similarly, n /= 2 is short for n = n / 2.

Exercise

Modify this algorithm to use repeated cubing instead of squaring and repeat this analysis. Is this approach faster than repeated squaring?