Skip to main content

Fibonacci sequence

ยท 3 min read
Josh Cena

This is the first activity of Computerization algorithm team. We introduced the method to find the nnth term of the Fibonacci sequence, which mainly uses matrix exponentiation.

Problemโ€‹

The Fibonacci sequence:

Fn={0,n=01,n=1Fnโˆ’2+Fnโˆ’1,n>1F_{n}= \begin{cases} 0,&n=0\\ 1,&n=1\\ F_{n-2}+F_{n-1},&n>1 \end{cases}

Given nn, find Fnย modย 109+7F_{n}\text{ mod }10^9+7ใ€‚

Input constraintsMemory limitExecution time
0โ‰คnโ‰ค10190\le n\le 10^{19}64MB1.0s

Solutionโ€‹

The input size of 101910^{19} obviously prohibits any attempt to solve it with loops. Is there a better way than a simple O(n)\mathcal{O}(n)? It turns out that with matrix exponentiation, we can achieve O(logโกn)\mathcal{O}(\log n). We observe that:

(Fn+1Fn+2)=(Fn+1Fn+Fn+1)=(0111)(FnFn+1)\begin{pmatrix}F_{n+1}\\F_{n+2}\end{pmatrix}=\begin{pmatrix}F_{n+1}\\F_{n}+F_{n+1}\end{pmatrix}=\begin{pmatrix}0&1\\1&1\end{pmatrix}\begin{pmatrix}F_{n}\\F_{n+1}\end{pmatrix}

This step is applicable to all recursive sequences, so it should be easily reached for an experienced candidate. Generally, for Fn+2=aFn+bFn+1F_{n+2}=aF_{n}+bF_{n+1}, we have

(Fn+1Fn+2)=(Fn+1aFn+bFn+1)=(01ab)(FnFn+1)\begin{pmatrix}F_{n+1}\\F_{n+2}\end{pmatrix}=\begin{pmatrix}F_{n+1}\\aF_{n}+bF_{n+1}\end{pmatrix}=\begin{pmatrix}0&1\\a&b\end{pmatrix}\begin{pmatrix}F_{n}\\F_{n+1}\end{pmatrix}

From the recursive definition,

(Fn+mFn+m+1)=(0111)m(FnFn+1)\begin{pmatrix}F_{n+m}\\F_{n+m+1}\end{pmatrix}=\begin{pmatrix}0&1\\1&1\end{pmatrix}^m\begin{pmatrix}F_{n}\\F_{n+1}\end{pmatrix}

Substituting n=0n=0,

(FmFm+1)=(0111)m(F0F1)\begin{pmatrix}F_{m}\\F_{m+1}\end{pmatrix}=\begin{pmatrix}0&1\\1&1\end{pmatrix}^m\begin{pmatrix}F_0\\F_1\end{pmatrix}

Now the problem is transformed into finding the matrix raised to the mmth power. If m=20a0+21a1+22a2+โ€ฆm=2^0a_0+2^1a_1+2^2a_2+\dots (representation in binary), then

(0111)m=((0111)1)a0ร—((0111)2)a1ร—((0111)4)a2โ€ฆ\begin{pmatrix}0&1\\1&1\end{pmatrix}^m=\left(\begin{pmatrix}0&1\\1&1\end{pmatrix}^{1}\right)^{a_0}\times \left(\begin{pmatrix}0&1\\1&1\end{pmatrix}^{2}\right)^{a_1}\times \left(\begin{pmatrix}0&1\\1&1\end{pmatrix}^{4}\right)^{a_2}\dots

The 2k2^kth powers of the original matrix can, in fact, be preprocessed. When m<1019m<10^{19}, k<logโก21019<64k<\log_2 10^{19}<64, so we only need to store at most 63 matrices. In addition,

(0111)2k=(0111)2kโˆ’1ร—(0111)2kโˆ’1\begin{pmatrix}0&1\\1&1\end{pmatrix}^{2^k}=\begin{pmatrix}0&1\\1&1\end{pmatrix}^{2^{k-1}}\times \begin{pmatrix}0&1\\1&1\end{pmatrix}^{2^{k-1}}

which implies that the powers can be attained within O(logโกm)\mathcal{O}(\log m) time. This is the idea of fast matrix exponentiation: compute all 2k2^kth powers, and put those needed together.

Programโ€‹

Below is the C++ code, where the most intractable part is probably implementation of matrix multiplication:

#include <iostream>
#include <cmath>

using namespace std;

struct mat {
unsigned long long a[4];
mat operator *(mat o) {
mat t;
t.a[0] = (this->a[0] * o.a[0] + this->a[1] * o.a[2]) % 1000000007;
t.a[1] = (this->a[0] * o.a[1] + this->a[1] * o.a[3]) % 1000000007;
t.a[2] = (this->a[2] * o.a[0] + this->a[3] * o.a[2]) % 1000000007;
t.a[3] = (this->a[2] * o.a[1] + this->a[3] * o.a[3]) % 1000000007;
return t;
}
};

// Preprocessed matrices raised to the 2^k power
mat mat_pow[64];

int fib(unsigned long long k) {
// Temporary matrix; each time multiply it by some term in mat_pow
mat tmp;
tmp.a[0] = 1;
tmp.a[1] = 0;
tmp.a[2] = 0;
tmp.a[3] = 1;
for (int i = 0; i < 64; i++) {
// If a_i is 1
if (k & (1ull << i)) {
tmp = tmp * mat_pow[i];
}
}
return tmp.a[1];
}

int main() {
mat_pow[0].a[0] = 0;
mat_pow[0].a[1] = 1;
mat_pow[0].a[2] = 1;
mat_pow[0].a[3] = 1;
for (int i = 1; i < 64; i++) {
mat_pow[i] = mat_pow[i-1] * mat_pow[i-1];
}
unsigned long long n;
cin >> n;
cout << fib(n) << endl;
return 0;
}

The formulae for matrix multiplication are:

(a0a1a2a3)ร—(b0b1b2b3)=(a0b0+a1b2a0b1+a1b3a2b0+a3b2a2b1+a3b3)\begin{pmatrix}a_0&a_1\\a_2&a_3\end{pmatrix}\times\begin{pmatrix}b_0&b_1\\b_2&b_3\end{pmatrix}=\begin{pmatrix}a_0b_0+a_1b_2&a_0b_1+a_1b_3\\a_2b_0+a_3b_2&a_2b_1+a_3b_3\end{pmatrix}

(a0a1a2a3)ร—(b0b1)=(a0b0+a1b1a2b0+a3b1)\begin{pmatrix}a_0&a_1\\a_2&a_3\end{pmatrix}\times\begin{pmatrix}b_0\\b_1\end{pmatrix}=\begin{pmatrix}a_0b_0+a_1b_1\\a_2b_0+a_3b_1\end{pmatrix}