Fibonacci in O(1) space and matrix exponentiation

Why this matters

The previous lesson stored one answer for every subproblem in a table. That works, but the table grows with the input. For Fibonacci, computing F(n) with a table of n + 1 numbers uses memory in proportion to n. Dynamic programming, usually shortened to DP, often has this shape: a table of answers that only needs its last few rows. Noticing that a table can shrink is a useful habit in DP.

This lesson does two things. First, it removes the table and keeps only two variables. That is O(1) space, which means the memory used does not grow with n. Second, it shows a different method for very large n. Instead of computing every Fibonacci number up to n, it multiplies a small matrix with itself, which needs only about log n steps. A matrix here is a grid of numbers, and you will see exactly what that means below.

You need to know how to write a for loop, how to use std::uint64_t from <cstdint>, and what a function parameter is. The memoization lesson introduced the Fibonacci recurrence, so this lesson uses the same definition: F(0) is 0, F(1) is 1, and F(n) is F(n - 1) + F(n - 2) for n of 2 or more.

Two terms are used throughout. Space means the memory a program needs while it runs. Time means how many steps the program takes, measured as the input grows. O(1) space means the memory is a fixed amount no matter how big the input is. O(log n) time means the number of steps grows slowly: doubling n adds only one more step.

Step 1: keep only the last two numbers

Look at how the table is filled. Each new row needs the two rows directly above it, and nothing older. So an old row can be thrown away as soon as it is used.

Here is F(n) for small n, with the last two values kept at each step. The variable prev holds F(i) and curr holds F(i + 1).

After each step, the pair moves one place forward. The new prev is the old curr, and the new curr is the sum. The loop runs n times, so the time is O(n). Only prev, curr and next exist, so the space is O(1).

Here is the code. It is the first function of the sample program that the rest of the lesson builds on.

#include <cstdint>
 
// Two variables: each step keeps only the last two Fibonacci numbers.
std::uint64_t fib_two_vars(int n) {
    std::uint64_t prev = 0;  // F(0)
    std::uint64_t curr = 1;  // F(1)
    for (int i = 0; i < n; ++i) {
        std::uint64_t next = prev + curr;
        prev = curr;
        curr = next;
    }
    return prev;
}

The loop runs zero times for n = 0, so the function returns prev, which is F(0). The loop runs once for n = 1. It sets next to 1, moves prev to 1, and returns prev, which is F(1). Check the value for n = 5 in the table above: after five steps, prev is 5.

The idea is the same as the table, with memory cut down. Each answer still depends on the two answers before it, so this is still dynamic programming. Only the storage has changed.

What the saving buys

Memory is small for Fibonacci numbers either way. The saving matters more in other DP problems, where each row is a long array. A table of size n by m can become two rows of size m. Later lessons in this module use that trick several times. For now, the important point is that the rows you no longer need can be dropped.

Step 2: write the step as a matrix

The two-variable loop moves a pair of numbers one place forward each time. That move can be written as a matrix multiplication, and this is the key to the second method.

A matrix is a rectangular grid of numbers. A 2 by 2 matrix has two rows and two columns. Write the current pair as a column: the top entry is F(n) and the bottom entry is F(n - 1). Write the step as a fixed 2 by 2 matrix M:

M = | 1  1 |
    | 1  0 |

Multiply M by the column [F(n), F(n - 1)]. The rule for multiplying is: each entry of the result is a row of M times the column, multiplied pair by pair and then added.

top row:    1 * F(n) + 1 * F(n - 1) = F(n) + F(n - 1) = F(n + 1)
bottom row: 1 * F(n) + 0 * F(n - 1) = F(n)

So M times [F(n), F(n - 1)] gives [F(n + 1), F(n)]. The matrix takes the pair one step forward, exactly like the loop did. That is the whole idea. The loop applied this step n times. A matrix can be applied n times by raising it to the power n.

Powers of M by hand

Multiplying M by itself gives M squared. Use the same rule: each entry is a row of the first matrix times a column of the second, pair by pair, then added.

M squared = | 1  1 | * | 1  1 | = | 1*1 + 1*1   1*1 + 1*0 | = | 2  1 |
            | 1  0 |   | 1  0 |   | 1*1 + 0*1   1*1 + 0*0 |   | 1  1 |

Multiply once more to get M cubed:

M cubed = M squared * M = | 2  1 | * | 1  1 | = | 2*1 + 1*1   2*1 + 1*0 | = | 3  2 |
                          | 1  1 |   | 1  0 |   | 1*1 + 1*1   1*1 + 1*0 |   | 2  1 |

Look at the entries. M to the power 3 is [[3, 2], [2, 1]]. The numbers 3, 2, 2, 1 are Fibonacci numbers, and the pattern holds in general:

M to the power n = | F(n + 1)   F(n)     |
                   | F(n)       F(n - 1) |

This pattern is easy to check for n = 1. The formula gives [[F(2), F(1)], [F(1), F(0)]], which is [[1, 1], [1, 0]]. That is M itself. For n = 3 it gives [[F(4), F(3)], [F(3), F(2)]], which is [[3, 2], [2, 1]]. The values match the hand calculation.

So the top-right entry of M to the power n is F(n). That is the entry the code will read.

Why the power helps

To get F(n), the code needs M to the power n. The obvious way is to multiply M by itself n times. That takes n multiplications, which is no better than the loop. The next step makes it fast.

Step 3: raise to a power by squaring

The trick is to build the power from squares, the same way a binary number is built from powers of two. Every whole number can be written as a sum of powers of two. The number 13 is 8 plus 4 plus 1, and in binary it is 1101.

So M to the power 13 can be written as:

M to the power 13 = M to the power 8 * M to the power 4 * M to the power 1

Each of these pieces comes from squaring the one before. Start with M and square it to get M to the power 2. Square that to get M to the power 4. Square that to get M to the power 8. Then multiply the three pieces needed for 13 together. That is three squarings and two multiplications, which is far less than 12 multiplications.

This method is called exponentiation by squaring. The algorithm reads the bits of the exponent from the lowest end. Each bit is either 0 or 1. Whenever the bit is 1, the current square is multiplied into the answer. Each bit also means one more squaring of the base.

Here is the process for exponent 13, written as a table. The bit column is the lowest bit of the exponent at that point. The base column is the power of M the loop is holding.

Read the rows as the loop runs. Each row halves the exponent by shifting it one bit to the right. The bit that falls off the end decides whether the current base is multiplied into the result. The base is squared on every row, so after k rows it is M to the power 2 to the power k.

The loop stops when the exponent reaches 0. The number of rows equals the number of bits in the exponent. The exponent n has about log2 n bits, so the loop runs about log2 n times. For n near one million, that is about 20 rows. For n near one quintillion (10 to the power 18), it is about 60 rows. Doubling n adds one more row.

The loop works on the exponent by shifting bits. You do not need to understand the binary arithmetic in depth. You only need to know that exp & 1 checks the lowest bit and exp >>= 1 removes it.

Step 4: the full program

Now put the pieces together. The matrix is a std::array of two std::array rows. Writing Matrix q{{{1, 1}, {1, 0}}}; gives q the values of M. The braces follow the matrix shape: an outer pair for the two rows, and an inner pair for each row's two numbers.

The next function multiplies two matrices. Step 5 explains the mod argument, which is the number the remainder is taken by. The power function uses the loop from Step 3. It starts with the identity matrix, the matrix that leaves other matrices unchanged, which is the matrix version of the number 1.

#include <array>
#include <cstdint>
#include <iostream>
 
// Two variables: each step keeps only the last two Fibonacci numbers.
std::uint64_t fib_two_vars(int n) {
    std::uint64_t prev = 0;  // F(0)
    std::uint64_t curr = 1;  // F(1)
    for (int i = 0; i < n; ++i) {
        std::uint64_t next = prev + curr;
        prev = curr;
        curr = next;
    }
    return prev;
}
 
using Matrix = std::array<std::array<std::uint64_t, 2>, 2>;
 
constexpr std::uint64_t MOD = 1'000'000'007ULL;
 
Matrix multiply(const Matrix& a, const Matrix& b, std::uint64_t mod) {
    Matrix c{};
    for (std::size_t i = 0; i < 2; ++i) {
        for (std::size_t j = 0; j < 2; ++j) {
            for (std::size_t k = 0; k < 2; ++k) {
                c[i][j] = (c[i][j] + a[i][k] * b[k][j]) % mod;
            }
        }
    }
    return c;
}
 
Matrix power(Matrix base, std::uint64_t exp, std::uint64_t mod) {
    Matrix result{{{1, 0}, {0, 1}}};  // identity matrix
    while (exp > 0) {
        if (exp & 1) {
            result = multiply(result, base, mod);
        }
        base = multiply(base, base, mod);
        exp >>= 1;
    }
    return result;
}
 
// F(n) mod `mod`, computed in O(log n) matrix multiplications.
std::uint64_t fib_matrix(std::uint64_t n, std::uint64_t mod) {
    Matrix q{{{1, 1}, {1, 0}}};
    return power(q, n, mod)[0][1];
}
 
int main() {
    std::cout << "n  F(n) two variables\n";
    for (int n = 0; n <= 10; ++n) {
        std::cout << n << "  " << fib_two_vars(n) << '\n';
    }
 
    std::cout << "\nmatrix power matches two variables for n = 0..90 (mod 1e9+7): ";
    bool all_match = true;
    for (int n = 0; n <= 90; ++n) {
        if (fib_matrix(n, MOD) != fib_two_vars(n) % MOD) {
            all_match = false;
            std::cout << "\nmismatch at n = " << n << '\n';
        }
    }
    std::cout << (all_match ? "yes" : "no") << '\n';
 
    std::cout << "\nF(10^18) mod 1e9+7 = " << fib_matrix(1'000'000'000'000'000'000ULL, MOD) << '\n';
    std::cout << "F(13) via matrix = " << fib_matrix(13, MOD) << '\n';
    return 0;
}

Compile and run it with the warnings turned on:

g++ -std=c++20 -Wall -Wextra -O2 -o fib fib.cpp
./fib

The real output is:

n  F(n) two variables
0  0
1  1
2  1
3  2
4  3
5  5
6  8
7  13
8  21
9  34
10  55
 
matrix power matches two variables for n = 0..90 (mod 1e9+7): yes
 
F(10^18) mod 1e9+7 = 209783453
F(13) via matrix = 233

The check on the line matrix power matches two variables compares the two methods for every n from 0 to 90. The line says yes, so both give the same answer for all of them. The last two lines show the matrix method on large and small inputs. Its value for F(13) is 233, which is F(13) in this numbering, where F(0) is 0.

The value for F(10^18) was computed by a second method in Python, which uses the same fast-doubling idea but not matrices. It gives 209783453 as well, so the two methods agree on that large input too. Python's integers can hold arbitrary sizes, so the Python check also confirms the 90 and 100 values below.

Nobody can check F(10^18) by hand. The check that counts is that the two-variable method agrees with the matrix method on every small input, and an independent method agrees on the huge one. Test a fast method against a simple one whenever you can.

Step 5: keep the numbers small with a modulus

The matrix method has a problem. Fibonacci numbers grow very fast. F(93) is the largest one that fits in a 64-bit unsigned integer, which holds values up to about 1.8 times 10 to the power 19. Beyond that, the value wraps around. This is called overflow. The program does not stop or print an error. It keeps going with a wrong number.

Here is the same method with the % mod step removed from multiply. The output for n = 90 is exact. For n = 100, it is not:

exact matrix power, n = 90: 2880067194370816120
exact matrix power, n = 100 (overflows 64 bits): 3736710778780434371

The true value of F(100) is 354224848179261915075. The program printed 3736710778780434371 instead, which is the true value with the 64-bit wrap applied. The output looks like a normal number, so nothing warns you. Overflow errors are quiet, and that is why they are dangerous.

Many problems do not want the exact Fibonacci number. They want the answer modulo a prime, such as 1'000'000'007. A modulus is the remainder after division. Writing x % mod gives the remainder of x divided by mod. When every step takes the remainder, the numbers never grow past the modulus, and the answer is still correct modulo that number. The multiply function does this with % mod after each addition.

The value 1,000,000,007 is a common choice for this reason. Each entry of a matrix stays below it. The product of two such numbers is below about 1.0 times 10 to the power 18, which fits in a 64-bit unsigned integer. If the modulus were much bigger, the products could overflow inside multiply.

Why the modulo answer is still correct

The rule that makes the modulo safe is that (a + b) % m and ((a % m) + (b % m)) % m give the same answer. The same holds for multiplication. So taking the remainder early, after every step, gives the same final remainder as taking it once at the end. The program never needs the huge intermediate values, so the result is correct modulo m.

Common mistakes

Forgetting the modulus. A program that multiplies big matrices without a modulus gives wrong answers from some point on, as the n = 100 output shows. The fix is to take the remainder inside multiply, after each product. Test with n near 100 to see the wrap.

Reading the wrong entry. The formula says the top-right entry of M to the power n is F(n). A common slip is to read [0][0] instead, which is F(n + 1). The test fib_matrix(13, MOD) == 233 catches this at once. Always check n = 0, n = 1 and n = 2 by hand when you write a matrix method.

Mixing up the identity matrix. The starting result must be the identity matrix [[1, 0], [0, 1]]. If you start with [[1, 1], [1, 0]], the answer is M to the power n + 1, which gives an off-by-one error in every answer. Start the result at the identity, and check n = 0, where the answer must be 0.

Using a signed type with a shift. The loop uses exp >>= 1. Shifting a negative number has its own rules, and a negative exponent makes no sense here. Use an unsigned type such as std::uint64_t for the exponent. The compiler warns when a signed and an unsigned value are compared. This is the message you get from -Wall when an int is compared with a std::uint64_t:

mistake.cpp:9:11: warning: comparison of integer expressions of different signedness: ‘int’ and ‘uint64_t’ {aka ‘long unsigned int’} [-Wsign-compare]
    9 |     if (n < count) {
      |         ~~^~~~~~~

Treat this warning as a bug to fix, not noise. Pick one type for the values and the counts, and convert once at the boundary.

Calling a function that lives in another file. Writing the helper in one file and calling it from another gives this error, which the compiler prints when the name is not visible:

mistake.cpp:4:18: error: ‘fib_two_vars’ was not declared in this scope
    4 |     std::cout << fib_two_vars(n) << '\n';
      |                  ^~~~~~~~~~~~

The fix is to declare the function in a header and include it, or to put both the function and its caller in the same file, as the sample program does.

Assuming matrix power is always faster. For a small n, such as 20, the loop does about 20 additions and the matrix method does several matrix multiplications, each with eight multiplications inside. The matrix method pays off only when n is large. Pick the method by the size of the input.

Practice

Exercise 1. Use the two-variable loop by hand to find F(20). Write the pair prev and curr after each step, or use the table pattern from Step 1.

Show the answer

F(20) is 6765. The pairs after each step are (0, 1), (1, 1), (1, 2), (2, 3), and so on. The sequence continues 0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597, 2584, 4181, 6765. Counting from F(0), the entry at index 20 is 6765.

Exercise 2. Use the pattern from Step 2 to find M to the power 4 without calculating it step by step from M. Then check it by multiplying M squared by itself.

Show the answer

The pattern says M to the power 4 is [[F(5), F(4)], [F(4), F(3)]], which is [[5, 3], [3, 2]]. To check, square M squared, which is [[2, 1], [1, 1]]. Each entry is a row times a column: the top-left is 2 times 2 plus 1 times 1, which is 5. The top-right is 2 times 1 plus 1 times 1, which is 3. The bottom-right is 1 times 1 plus 1 times 1, which is 2. The result matches.

Exercise 3. How many matrix multiplications does power do when n is 2 to the power 20, which is 1,048,576? Count the squarings, and count the result multiplications.

Show the answer

The exponent 2 to the power 20 has 21 binary digits: a single 1 followed by 20 zeros. The loop does one squaring for each bit, so 21 squarings in total. Only one bit is 1, so the result is multiplied only once. The total is 22 multiplications, which is far fewer than the 1,048,576 needed by multiplying M by itself in a loop. Note that the last squaring of the base is wasted in this code. A stricter version stops squaring when the exponent reaches 0, which saves one multiplication.

Exercise 4. Write a function fib_last_digit(n) that returns the last decimal digit of F(n), for example by passing 10 as the modulus to fib_matrix. Find the last digit of F(100) by running it.

Show the answer

Calling fib_matrix(100, 10) returns 5. The true value F(100) is 354224848179261915075, which ends in 5. The function is just return fib_matrix(n, 10);. This works because the remainder of a number divided by 10 is its last digit, and the modulo rule from Step 5 keeps the last digit correct through every multiplication.

Exercise 5. The two-variable method uses two variables for any recurrence where each value depends on the previous two. Write a function tribonacci_two_vars for the sequence T(0) = 0, T(1) = 0, T(2) = 1, and T(n) = T(n - 1) + T(n - 2) + T(n - 3). Which number of variables do you need, and why?

Show the answer

Three variables are needed, one for each of the three previous values. Each step computes the new value from a, b and c, then shifts them: a becomes b, b becomes c, and c becomes the new value. In general, a recurrence that looks back k places needs k variables. The space is still O(1) because k is fixed by the recurrence, not by n.

Summary

• Fibonacci needs only the last two values, so a loop with two variables uses O(1) space and O(n) time.
• The step from one pair to the next is a 2 by 2 matrix, and its power n holds F(n) in the top-right entry.
• Exponentiation by squaring computes a matrix power in about log n multiplications, using the bits of the exponent.
• Take the remainder modulo a prime after every multiplication, or the numbers overflow and the output is wrong without any warning.
• Check a fast method against a simple one on small inputs before trusting it on large ones.