【矩阵快速幂】洛谷 P1962 斐波那契数列

题目

https://www.luogu.com.cn/problem/P1962

题解

矩阵快速幂是一种利用快速幂思想来计算矩阵幂次的高效算法，它能将矩阵幂运算的时间复杂度从 \(O(n^3k)\) 降低到 \(O(n^3logk)\)，其中 \(n\) 是矩阵维度，\(k\) 是指数。

斐波那契数列（Fibonacci Sequence），又称为黄金分割数列，是一个经典的整数序列，其定义如下：

\[F_i = \begin{cases} 1, & i = 1 \text{ 或 } i = 2 \\ F_{i-1} + F_{i-2}, & i > 2 \end{cases} \]

斐波那契数列具有明显重复的递推关系，因此可以考虑使用线性代数的知识，将递推方程表示为矩阵乘积的形式：

\[\begin{bmatrix} F_{n-1} & F_{n} \end{bmatrix} = \begin{bmatrix} F_{n-2} & F_{n-1} \end{bmatrix} \cdot M = \begin{bmatrix} F_{n-2} & F_{n-1} \end{bmatrix} \cdot \begin{bmatrix} a & b \\ c & d \end{bmatrix} = \begin{bmatrix} a \cdot F_{n-2} + c \cdot F_{n-1} & b \cdot F_{n-2} + d \cdot F_{n-1} \end{bmatrix} \]

那么，得到以下关系式

\[\begin{cases} F_{n-1} = a \cdot F_{n-2} + c \cdot F_{n-1} \\ F_{n} = b \cdot F_{n-2} + d \cdot F_{n-1} \end{cases} \]

解方程，可得：

\[\begin{cases} a = 0, \\ b = c = d = 1. \end{cases} \]

亦即

\[M = \begin{bmatrix} 0 & 1 \\ 1 & 1 \end{bmatrix} \]

于是可得求第 \(k\) 阶斐波那契数列的矩阵乘积为：

\[\begin{bmatrix} F_{k-1} & F_{k} \end{bmatrix} = \begin{bmatrix} F_1 & F_2 \end{bmatrix} \cdot M^{k-2} = \begin{bmatrix} 1 & 1 \end{bmatrix} \cdot \begin{bmatrix} 0 & 1 \\ 1 & 1 \end{bmatrix}^{k-2} \]

综上，可以使用矩阵快速幂加速 \(M^{k-2}\) 的运算，从而达到快速求出斐波那契数列的第 \(k\) 项的值。

参考代码

#include<bits/stdc++.h>
using namespace std;
using ll = long long;

constexpr int MOD = 1e9 + 7;

template<typename T>
class MatrixQPow {
private:
    using MATRIX_TYPE = vector<vector<T>>;

    MATRIX_TYPE multiple(const MATRIX_TYPE &matrix1, const MATRIX_TYPE &matrix2, const ll mod = LONG_LONG_MAX) {
        int n1 = matrix1.size(), m1 = matrix1[0].size();
        int n2 = matrix2.size(), m2 = matrix2[0].size();
        MATRIX_TYPE res(n1, vector<T>(m2, 0));// 默认视 m1 == n2 成立，可生成 n1 × m2 的矩阵
        for (int i = 0; i < n1; ++i) {
            for (int j = 0; j < m2; ++j) {
                for (int k = 0; k < m1; ++k) {
                    res[i][j] = (res[i][j] + matrix1[i][k] * matrix2[k][j]) % mod;
                }
            }
        }
        return res;
    }

public:
    /*矩阵快速幂*/
    MATRIX_TYPE qpow(MATRIX_TYPE &base, ll exponent, const ll mod = LONG_LONG_MAX) {
        int n = base.size();
        MATRIX_TYPE res(n, vector<T>(n, 0));
        for (int i = 0; i < n; ++ i) res[i][i] = 1;// 生成单位矩阵
        while (exponent > 0) {
            if (exponent & 1) res = multiple(res, base, mod);
            base = multiple(base, base, mod);
            exponent >>= 1;
        }
        return res;
    }

    /*Fibonacci 数列求第 n 项*/
    T getNthFibonacci(ll n, const ll mod = LONG_LONG_MAX) {
        if (n <= 0) return 0;
        if (n == 1 || n == 2) return 1;

        // 转移矩阵 M = [[0,1],[1,1]]
        MATRIX_TYPE m = {
                {0, 1},
                {1, 1}
        };
        // 初始向量 [F1, F2] = [1, 1]
        MATRIX_TYPE w = {{1, 1}};  // 明确使用1x2矩阵

        // 计算 m^(n-2)
        MATRIX_TYPE power_mat = qpow(m, n - 2, mod);
        // [F1, F2] * m^(n-2) = [F(n-1), F(n)]
        MATRIX_TYPE fib = multiple(w, power_mat, mod);
        return fib[0][1];  // 取F(n)
    }
};

int main() {
    ios::sync_with_stdio(false);cin.tie(nullptr);cout.tie(nullptr);
    ll n;
    cin >> n;
    MatrixQPow<ll> matrixQPow;
    cout << matrixQPow.getNthFibonacci(n, MOD) << '\n';
    return 0;
}