# 4 cách khác nhau để tính số Fibonacci/truy hồi tuyến tính

Cho bài toán sau:

> Ta định nghĩa $F(i)$ là số Fibonacci thứ $i$ theo công thức truy hồi sau:
> $$F(i) = \begin{cases} 1 & (i \le 1) \\ F(i - 1) + F(i - 2) &(i \ge 2) \end{cases}.$$
> Tính $F(n)$ modulo $998244353$.

Mình sẽ giới thiệu 4 cách khác nhau để giải bài toán này, từ cách đơn giản nhất tới cách nâng cao nhất.

## Đệ quy

Một trong những cách đơn giản nhất là sử dụng đệ quy để trực tiếp cài đặt công thức truy hồi trên:

int F(int i) {
    return i <= 1 ? 1 : (F(i - 1) + F(i - 2)) % MOD;
}

Tuy nhiên, khi chạy thử code này, các bạn sẽ thấy code này không chạy được với $n \ge 50$. Lý do là vì cây đệ quy của cách cài đặt này cực kì to (nếu các bạn để ý, cây đệ quy này có độ lớn ít nhất bằng chính số $F(n)$ các bạn đang tính). Cũng chính vì lý do này mà độ phức tạp của cách làm này là $O(F(n)) = O(1.618^n)$.

Nhận xét rằng code trên có thể sẽ gọi rất nhiều lần $F(k)$ (với $k \le n$) làm code bị chậm. Tuy nhiên, đáp án của $F(k)$ không bao giờ thay đổi dù ta có gọi $F(k)$ bao nhiêu lần. Vì thế, _nếu_ ta lưu đáp án của $F(k)$ rồi trả về đáp án này mỗi lần gọi sau, ta có thể cải tiến thời gian rất nhiều.

## Quy hoạch động (và đệ quy có nhớ)

Nhận xét trên cho ta cách cài đặt sau.

vector<int> fib(n + 1, -1);
int F(int i) {
    if (fib[i] != -1) {
        return fib[i];
    } else {
        fib[i] = i <= 1 ? 1 : (F(i - 1) + F(i - 2)) % MOD;
        return fib[i];
    }
}

Cách cài đặt này cho ta độ phức tạp là $O(n)$, và với cách làm này ta có thể tính đến số Fibonacci thứ $10^8$ trong khoảng 1 giây! Cách cài đặt này được gọi là đệ quy có nhớ (hoặc là quy hoạch động top-down), và đây là một cách để cải tiến các bài toán khác có nhận xét tương tự. Tuy nhiên, đây không phải là cách phổ biến nhất để cài đặt quy hoạch động; ta có thể dùng vòng for để tính đáp án như sau.

vector<int> fib(n + 1);
fib[0] = fib[1] = 1;
for (int i = 2; i <= n; i++) {
    fib[i] = (fib[i - 1] + fib[i - 2]) % MOD;
}

Tuy nhiên, ta không thể dùng cách này để giải bài toán với $n$ lớn hơn (ví dụ như $n \le 10^{18}$).

## Nhân ma trận

Để ý rằng công thức quy hoạch động trên "đồng nhất": cách tính $F(i)$ dựa trên $F(i - 1)$ và $F(i - 2)$ giống hệt cách tính $F(k)$ dựa trên $F(k - 1)$ và $F(k - 2)$. Tổng quát hơn, giả sử ta lưu hai số Fibonacci được tính cuối cùng thành vectơ $\begin{bmatrix} F(i) \\ F(i - 1) \end{bmatrix}$, ta có thể tính vectơ $\begin{bmatrix} F(i + 1) \\ F(i) \end{bmatrix}$ như sau:
$$\begin{bmatrix} F(i + 1) \\ F(i) \end{bmatrix} = \begin{bmatrix} F(i) + F(i - 1) \\ F(i) \end{bmatrix} = \begin{bmatrix}1 & 1 \\ 1 & 0\end{bmatrix} \begin{bmatrix} F(i) \\ F(i - 1) \end{bmatrix}.$$

Vì thế, ta có thể chứng minh bằng đệ quy dựa trên nhận xét trên rằng $$\begin{bmatrix} F(n) \\ F(n - 1) \end{bmatrix} = \begin{bmatrix}1 & 1 \\ 1 & 0\end{bmatrix}^{n - 1} \begin{bmatrix} F(1) \\ F(0) \end{bmatrix} = \begin{bmatrix}1 & 1 \\ 1 & 0\end{bmatrix}^{n - 1} \begin{bmatrix} 1 \\ 1 \end{bmatrix}.$$

Việc tính nhanh $A^k$ với số $k$ bất kì có thể được thực hiện bằng thuật toán tính lũy thừa nhanh như sau:
$$A^k = \begin{cases} (A^{\lfloor k/2 \rfloor})^2 & \textrm{($k$ chẵn)} \\ (A^{\lfloor k/2 \rfloor})^2 \cdot A & \textrm{($k$ lẻ)} \end{cases}.$$

Và vì thế, ta có thuật toán tính $F(n)$ với độ phức tạp là $O(\log n)$. Kĩ thuật nhân ma trận này được sử dụng cực kì nhiều cho các bài toán quy hoạch động có công thức truy hồi đồng nhất, với độ phức tạp chung là $O(|A|^3 \log n)$, với $|A|$ là kích cỡ của công thức truy hồi và $n$ là đáp án cần tính.

## Bí thuật hàm sinh và thuật toán Bostan-Mori
Ta quay lại bài toán ban đầu được khái quát thành bài toán truy hồi tuyến tính như sau:
> Cho mảng $c$ có kích cỡ là $k$. Ta định nghĩa $F(i)$ là số $k$-Fibonacci thứ $i$ theo công thức truy hồi sau:
> $$F(i) = \begin{cases} 1 & (0 \le i < k) \\ \sum_{i = j}^k c_j \cdot F(i - j) &(i \ge k) \end{cases}.$$
> Tính $F(n)$ modulo $998244353$.

Với $k \le 300$, ta hoàn toàn có thể giải bài toán trên với độ phức tạp là $O(k^3 \log n)$ sử dụng nhân ma trận. Tuy nhiên, với $k \le 10^5$, ta cần phải tìm một thuật toán khác. Ta có thể dùng hàm sinh để giải bài toán này. Với các bạn chưa biết, với một dãy số $a$ bất kì, hàm sinh $g_a$ của $a$ được định nghĩa là chuỗi lũy thừa (power series) sau: $g_a(x) = \sum_{i = 0}^{\infty} a_i x^i$.

Với dãy $c$, ta đặt $Q(x) = 1 - c_{1} x - c_{2} x^2 - \dots - c_{k} x^k = 1 - \sum_{i = 1}^{k} c_{j} x^{j}$. Nhận xét rằng nếu ta nhân hàm sinh của $F$ cho $Q$, mọi lũy thừa $x^i$ với $i \ge k$ của $g_F(x) \cdot Q(x)$ đều có hệ số bằng $0$:
$$[x^i]g_F(x) \cdot Q(x) = \sum_{j=0}^k [x^{i-j}]g_F(x) \cdot [x^{j}]Q(x) = F(i) - \sum_{j=1}^k F(i - j) \cdot c_j = 0.$$

Vì thế, $P(x) := Q(x) \cdot g_F(x)$ là một đa thức với bậc bé hơn $k$. Ngoài ra, ta hoàn toàn có thể biết được hệ số của $x^0$ tới $x^{k - 1}$ của $P(x)$ bằng cách nhân đa thức $\sum_{i = 0}^{k - 1} F(i) x^i = \sum_{i = 0}^{k - 1} x^i$ cho $Q(x)$ (bởi vì trong $g_F(x)$, chỉ có các hệ số từ $x^0$ tới $x^{k - 1}$ đóng góp vào $k$ hệ số đầu của $P(x)$). Vì thế, ta có thể tính được chính xác $P(x)$ và $Q(x)$ bằng nhân đa thức nhanh trong $O(k \log k)$. Bởi vì $g_F(x) = \frac{P(x)}{Q(x)}$, bài toán được chuyển về như sau:

> Cho đa thức $P(x)$ và $Q(x)$ với bậc tối đa là $k$. Tính hệ số của $x^n$ của $\frac{P(x)}{Q(x)}$.

Ý tưởng cơ bản là như sau: nhận xét rằng $Q(x) Q(-x)$ là hàm chẵn, vì thế nên mọi hệ số ở các lũy thừa chẵn của $Q(x) Q(-x)$ là bằng $0$. Chính vì thế, nếu ta đặt $Q(x) Q(-x) = V(x^2)$ với đa thức $V$ nào đó, thì
$$\frac{P(x)}{Q(x)} = \frac{P(x) Q(-x)}{Q(x) Q(-x)} = \frac{U_e(x^2)}{V(x^2)} + x \cdot \frac{U_o(x^2)}{V(x^2)},$$
với $U_e$ và $U_o$ lần lượt là đa thức gồm các hệ số ở lũy thừa chẵn và lẻ của $U(x) = P(x) Q(-x)$. Từ công thức trên, nhận xét rằng
$$[x^n]\frac{P(x)}{Q(x)} = \begin{cases} [x^n]\frac{U_e(x^2)}{V(x^2)} = [x^{\lfloor n/2 \rfloor }]\frac{U_e(x)}{V(x)} & \textrm{($n$ chẵn)} \\ [x^n]\frac{x U_o(x^2)}{V(x^2)} = [x^{\lfloor n/2 \rfloor }]\frac{U_o(x)}{V(x)} & \textrm{($n$ lẻ)} \end{cases}.$$

Vì thế, ta có thể đưa bài toán "tìm hệ số của $x^n$" thành một bài toán tìm hệ số khác với $n$ giảm đi một nửa sau hai lần nhân đa thức nhanh. Với base case là $[x^0]\frac{P(x)}{Q(x)} = \frac{P(0)}{Q(0)}$, ta có thể giải bài toán trên với độ phức tạp là $O(k \log k \log n)$.

Thuật toán trên được gọi là thuật toán Bostan-Mori ở cộng đồng Nhật Bản; các bạn có thể tìm paper của thuật toán này [tại đây](https://arxiv.org/pdf/2008.08822.pdf). Tuy nhiên, kĩ thuật tính hệ số trên không chỉ áp dụng cho mỗi bài toán tính truy hồi tuyến tính nhanh mà còn có thể áp dụng cho nhiều bài toán về hàm sinh khác.

## Một số bài tập
Nhân ma trận: [Problem B - NAQ 2022](https://naq22.kattis.com/contests/naq22/problems/birthdaygift).
Bostan-Mori: [Problem T - AtCoder Typical DP Contest](https://atcoder.jp/contests/tdpc/tasks/tdpc_fibonacci).