# Fibonacci sequence

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

## Problem​

The Fibonacci sequence:
$F_{n}= \begin{cases} 0,&n=0\\ 1,&n=1\\ F_{n-2}+F_{n-1},&n>1 \end{cases}$

Given $n$, find $F_{n}\text{ mod }10^9+7$

Input constraintsMemory limitExecution time
$0\le n\le 10^{19}$64MB1.0s

## Solution​

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

$\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 $F_{n+2}=aF_{n}+bF_{n+1}$, we have

$\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,

$\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=0$,

$\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 $m$th power. If $m=2^0a_0+2^1a_1+2^2a_2+\dots$ (representation in binary), then

$\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 $2^k$th powers of the original matrix can, in fact, be preprocessed. When $m<10^{19}$, $k<\log_2 10^{19}<64$, so we only need to store at most 63 matrices. In addition,

$\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 $\mathcal{O}(\log m)$ time. This is the idea of fast matrix exponentiation: compute all $2^k$th 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 powermat 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:

$\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}$

$\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}$
{"extraction_info": {"found_math": true, "script_math_tex": 0, "script_math_asciimath": 0, "math_annotations": 25, "math_alttext": 0, "mathml": 0, "mathjax_tag": 0, "mathjax_inline_tex": 0, "mathjax_display_tex": 0, "mathjax_asciimath": 1, "img_math": 0, "codecogs_latex": 0, "wp_latex": 0, "mimetex.cgi": 0, "/images/math/codecogs": 0, "mathtex.cgi": 0, "katex": 0, "math-container": 0, "wp-katex-eq": 0, "align": 0, "equation": 0, "x-ck12": 0, "texerror": 0, "math_score": 0.7334580421447754, "perplexity": 2863.598544482607}, "config": {"markdown_headings": true, "markdown_code": true, "boilerplate_config": {"ratio_threshold": 0.18, "absolute_threshold": 10, "end_threshold": 15, "enable": true}, "remove_buttons": true, "remove_image_figures": true, "remove_link_clusters": true, "table_config": {"min_rows": 2, "min_cols": 3, "format": "plain"}, "remove_chinese": true, "remove_edit_buttons": true, "extract_latex": true}, "warc_path": "s3://commoncrawl/crawl-data/CC-MAIN-2021-43/segments/1634323585522.78/warc/CC-MAIN-20211022212051-20211023002051-00181.warc.gz"}