矩阵快速幂学习笔记

矩阵快速幂学习笔记

首先我们得先知道什么是快速幂

这里简要提一下，就是求 \(a^n\) 次方时，我们可以把 \(n\) 拆分成二进制(假设为 \(.....0110\))，然后就相当于求 \(a^x+a^y......a^4+a^2\) 。

然后矩阵快速幂就是对一个矩阵做快速幂，好像是废话

其实矩阵快速幂的实现非常简单，就是将快速幂中的乘法变为矩阵乘法即可，这里重点讲用矩阵快速幂优化递推。

模板：

vector<vector<LL>> mul(vector<vector<LL>> a, vector<vector<LL>> b, LL mod)
{
    vector<vector<LL>> res(a.size(), vector<LL>(b[0].size(), 0));
    int n = a.size(), m = b[0].size(), x = b.size();
    for (int i = 0; i < n; i++)
        for (int j = 0; j < m; j++)
            for (int k = 0; k < x; k++)
                res[i][j] = (res[i][j] + a[i][k] * b[k][j] % mod) % mod;
    return res;
}

vector<vector<LL>> qmi(vector<vector<LL>> a, LL k, LL mod)
{
    vector<vector<LL>> I = {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}};
    while (k)
    {
        if (k & 1)
            I = mul(I, a, mod);
        a = mul(a, a, mod);
        k >>= 1;
    }
    return I;
}

直接看例题：

P1962 斐波那契数列

易得答案矩阵为 \(\begin{bmatrix} f_n\\ f_{n-1} \end{bmatrix}\)

然后 \(f_n\) 是从 \(\begin{bmatrix} f_{n-1} \\ f_{n-2} \end{bmatrix}\) 根据矩阵乘法的性质，递推矩阵为 \(\begin{bmatrix} 1 & 1\\ 1 & 0 \end{bmatrix}\), 因为 \(f_n=f_{n-1}+f_{n-2}\)

那么现在我们知道 \(\begin{bmatrix} f_2 \\ f_1 \end{bmatrix}\) 就要将乘法进行 \(n-2\) 次，用矩阵快速幂优化。

所以递推式子就为：

\[\begin{bmatrix} 1 & 1\\ 1 & 0 \end{bmatrix}^{n-2}*\begin{bmatrix} f_2 \\ f_1 \end{bmatrix} \]

#include <algorithm>
#include <cstring>
#include <iostream>
#include <vector>
using namespace std;
typedef long long LL;
typedef pair<int, int> PII;
const LL mod = 1e9 + 7;
inline LL read()
{
    LL x = 0, f = 1;
    char s = getchar();
    while (s < '0' || s > '9')
    {
        if (s == '-')
            f = -f;
        s = getchar();
    }
    while (s >= '0' && s <= '9')
    {
        x = (x << 3) + (x << 1) + (s ^ 48);
        s = getchar();
    }
    return x * f;
}
//矩阵乘法
vector<vector<LL>> mul(vector<vector<LL>> A, vector<vector<LL>> B, LL MOD)
{
    vector<vector<LL>> res(A.size(), vector<LL>(B[0].size(), 0));
    int n = A.size(), m = B[0].size(), x = B.size();
    for (int i = 0; i < n; i++)
        for (int j = 0; j < m; j++)
            for (int k = 0; k < x; k++)
                res[i][j] = (res[i][j] + A[i][k] * B[k][j] % MOD) % MOD;
    return res;
}
//矩阵快速幂
vector<vector<LL>> qmi(vector<vector<LL>> A, LL k, LL MOD)
{
    //初始化单位矩阵
    int n = A.size();
    vector<vector<LL>> I(n, vector<LL>(n, 0));
    for (int i = 0; i < n; i++)
        I[i][i] = 1;
    while (k)
    {
        if (k & 1)
            I = mul(I, A, MOD);
        A = mul(A, A, MOD);
        k >>= 1;
    }
    return I;
}

int main()
{
    LL n;
    n = read();
    if (n <= 2)
    {
        cout << '1';
        return 0;
    }
    vector<vector<LL>> A{{1, 1}, {1, 0}}; //递推矩阵
    vector<vector<LL>> f3{{1}, {1}};
    auto M = qmi(A, n - 2, mod);
    auto fn = mul(M, f3, mod);
    cout << fn[0][0] << endl;
    return 0;
}

P1939 【模板】矩阵加速（数列）

易得答案矩阵为 $\begin{bmatrix}
f_{n} \ f_{n-1} \ f_{n-2}
\end{bmatrix} $

然后 \(f_n\) 是从推出来的 \(\begin{bmatrix} f_{n-1} \\ f_{n-2} \\ f_{n-3} \end{bmatrix}\) 根据矩阵乘法的性质，递推矩阵为 \(\begin{bmatrix} 1 & 0 & 1\\ 1 & 0 & 0\\ 0 & 1 & 0\end{bmatrix}\), 因为 \(f_n=f_{n-1}+f_{n-3}\)

那么现在我们知道 \(\begin{bmatrix} f_3 \\ f_2 \\ f_1 \end{bmatrix}\) 就要将乘法进行 \(n-3\) 次，用矩阵快速幂优化。

所以递推式子就为：

\[\begin{bmatrix} 1 & 0 & 1\\ 1 & 0 & 0\\ 0 & 1 & 0\end{bmatrix}^{n-3}*\begin{bmatrix} f_3 \\ f_2 \\ f_1 \end{bmatrix} \]

#include <algorithm>
#include <cstring>
#include <iostream>
#include <vector>
using namespace std;
typedef long long LL;
typedef pair<int, int> PII;
const LL mod = 1e9 + 7;
inline int read()
{
    int x = 0, f = 1;
    char s = getchar();
    while (s < '0' || s > '9')
    {
        if (s == '-')
            f = -f;
        s = getchar();
    }
    while (s >= '0' && s <= '9')
    {
        x = (x << 3) + (x << 1) + (s ^ 48);
        s = getchar();
    }
    return x * f;
}

vector<vector<LL>> mul(vector<vector<LL>> a, vector<vector<LL>> b, LL mod)
{
    vector<vector<LL>> res(a.size(), vector<LL>(b[0].size(), 0));
    int n = a.size(), m = b[0].size(), x = b.size();
    for (int i = 0; i < n; i++)
        for (int j = 0; j < m; j++)
            for (int k = 0; k < x; k++)
                res[i][j] = (res[i][j] + a[i][k] * b[k][j] % mod) % mod;
    return res;
}

vector<vector<LL>> qmi(vector<vector<LL>> a, LL k, LL mod)
{
    vector<vector<LL>> I = {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}};
    while (k)
    {
        if (k & 1)
            I = mul(I, a, mod);
        a = mul(a, a, mod);
        k >>= 1;
    }
    return I;
}

int main()
{
    int q;
    cin >> q;
    while (q -- )
    {
        LL n = read();
        if (n <= 3)
        {
            cout << 1 << endl;
            continue;
        }
        vector<vector<LL>> a = {{1, 0, 1}, {1, 0, 0}, {0, 1, 0}};
        vector<vector<LL>> f3 = {{1}, {1}, {1}};
        vector<vector<LL>> M = qmi(a, n - 3, mod);
        auto ans = mul(M, f3, mod);
        cout << ans[0][0] << endl;
    }
    return 0;
}