2021 暑假训练学习矩阵快速幂

矩阵快速幂

矩阵快速幂的引入

给定一个大小为 \(n\) 的方阵 \(A\)，求 \(A^k\)。

\(1\leq n \leq 100,0\leq k \leq 10^{12},|A_{i,j}|\leq 1000\)

题目来源：P3390 【模板】矩阵快速幂

令 \(C=A*B\)，则 \(C_{i,j}=\sum\limits_{k=1}^{n}A_{i,k}B_{k,j}\) ，显然每次乘法的复杂度是 \(O(n^3)\) 的。

如果我们只是普通的进行递推乘法，那么复杂度是 \(O(n^3k)\) 的，显然无法承受。

不过，我们可以使用快速幂的方式来优化复杂度（即 \(A^{k}=(A^{\frac{k}{2}})^2\)），从而优化到 \(O(n^3 \log k)\) 。

#include<bits/stdc++.h>
using namespace std;
#define LL long long
const int N = 110;
const LL mod = 1e9 + 7;
struct Matrix {
    int n;
    LL a[N][N];
    void init(int _n, int opt) {
        memset(a, 0, sizeof(a));
            n = _n;
            if (opt == 1)
                for (int i = 1; i <= n; ++i)
                    a[i][i] = 1;
    }
    Matrix operator = (const Matrix &B) {
        this->n = B.n;
        for (int i = 1; i <= n; ++i)
            for (int j = 1; j <= n; ++j)
                this->a[i][j] = B.a[i][j];
        return *this;
    }
    Matrix operator * (const Matrix &B) const {
        Matrix res;
        res.init(n, 0);
        for (int i = 1; i <= n; ++i)
            for (int j = 1; j <= n; ++j)
                for (int k = 1; k <= n; ++k) {
                    res.a[i][j] += (this->a[i][k]) * B.a[k][j];
                    res.a[i][j] %= mod;
                }
        return res;
    }
    void output() {
        for (int i = 1; i <= n; ++i) {
            for (int j = 1; j < n; ++j)
                printf("%lld ", a[i][j]);
            printf("%lld\n", a[i][n]);
        }
    }
};
Matrix quickpow(Matrix P, LL n) {
    if (n == 1) return P;
    Matrix res;
    res.init(P.n, 1);
    while (n) {
        if (n & 1) res = res * P;
        n >>= 1;
        P = P * P;
    }
    return res;
}
int main()
{
    int n;
    LL k;
    cin >> n >> k;
    Matrix A;
    A.init(n, 0);
    for (int i = 1; i <= n; ++i)
        for (int j = 1; j <= n; ++j)
            cin >> A.a[i][j];
    Matrix B = quickpow(A, k);
    B.output();
    return 0;
}

应用：优化线性递推式

求斐波那契数列的第 \(n\) 项 \(f_n\)。（\(f_1=f_2=1,f_n=f_{n-1}+f_{n-2}(n \geq 3)\)）

\(1\leq n \leq 2^{63}\)

题目来源：P1962 斐波那契数列

递推的复杂度是 \(O(n)\)，现在不能接受。

所有递推式都可以表示成矩阵乘法的形式，例如

\[f_n= \begin{bmatrix} 1&1 \end{bmatrix} \begin{bmatrix} f_{n-1} \\ f_{n-2} \end{bmatrix} \]

不过这个显然没法进行矩阵快速幂，所以我们改一下，变成：

\[\begin{bmatrix} f_{n} \\ f_{n-1} \end{bmatrix} = \begin{bmatrix} 1&1 \\ 1&0 \end{bmatrix} \begin{bmatrix} f_{n-1} \\ f_{n-2} \end{bmatrix} \]

对于 \(n\geq 3\)，有

\[\begin{bmatrix} f_{n} \\ f_{n-1} \end{bmatrix} = \begin{bmatrix} 1&1 \\ 1&0 \end{bmatrix} \begin{bmatrix} f_{n-1} \\ f_{n-2} \end{bmatrix} = \begin{bmatrix} 1&1 \\ 1&0 \end{bmatrix} ^2 \begin{bmatrix} f_{n-2} \\ f_{n-3} \end{bmatrix} =\cdots= \begin{bmatrix} 1&1 \\ 1&0 \end{bmatrix} ^{n-2} \begin{bmatrix} f_{2} \\ f_{1} \end{bmatrix} \]

有初状态 \(f_1=f_2=1\)，所以我们直接矩阵快速幂求出后面的矩阵，和初状态矩阵乘一下即可得出答案。

#include<bits/stdc++.h>
using namespace std;
#define LL long long
const int N = 5;
const LL mod = 1e9 + 7;
struct Matrix {
    int n;
    LL a[N][N];
    void init(int _n, int opt) {
        memset(a, 0, sizeof(a));
            n = _n;
            if (opt == 1)
                for (int i = 1; i <= n; ++i)
                    a[i][i] = 1;
    }
    Matrix operator = (const Matrix &B) {
        this->n = B.n;
        for (int i = 1; i <= n; ++i)
            for (int j = 1; j <= n; ++j)
                this->a[i][j] = B.a[i][j];
        return *this;
    }
    Matrix operator * (const Matrix &B) const {
        Matrix res;
        res.init(n, 0);
        for (int i = 1; i <= n; ++i)
            for (int j = 1; j <= n; ++j)
                for (int k = 1; k <= n; ++k) {
                    res.a[i][j] += (this->a[i][k]) * B.a[k][j];
                    res.a[i][j] %= mod;
                }
        return res;
    }
};
Matrix quickpow(Matrix P, LL n) {
    if (n == 1) return P;
    Matrix res;
    res.init(P.n, 1);
    while (n) {
        if (n & 1) res = res * P;
        n >>= 1;
        P = P * P;
    }
    return res;
}
void output(const Matrix &ANS) {
    LL ans = ANS.a[1][1] + ANS.a[1][2] + ANS.a[1][3];
    cout << ans % mod << endl;
}
int main()
{
    //init
    Matrix A;
    A.init(2, 0);
    A.a[1][1] = A.a[1][2] = A.a[2][1] = 1;
    //
    LL n;
    cin >> n;
    if (n <= 2) {
        puts("1");
        return 0;
    }
    Matrix P = quickpow(A, n - 2);
    output(P);
    return 0;
}

不过不是所有递推式都这么简单，所以我们列出一些式子，来看看如何转换：

\[f_n=2f_{n-2}+f_{n-1}+1 \\ \begin{bmatrix} f_{n} \\ f_{n-1} \\ 1 \end{bmatrix} = \begin{bmatrix} 1 & 2 & 1 \\ 1 & 0 & 0 \\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} f_{n-1} \\ f_{n-2} \\ 1 \end{bmatrix} \]

下面这个很奇妙，因为引入了除了 \(f\) 外，新的可变数列（这时我们需要导入常数来维护它）：

\[f_n=7f_{n-1}+6f_{n-2}+5n+4*3^n \\ \begin{bmatrix} f_{n}\\f_{n-1}\\n\\3^n\\1 \end{bmatrix} = \begin{bmatrix} 7&6&5&12&5\\ 1&0&0&0&0\\ 0&0&1&0&1\\ 0&0&0&3&0\\ 0&0&0&0&1 \end{bmatrix} \begin{bmatrix} f_{n-1}\\f_{n-2}\\n-1\\3^{n-1}\\1 \end{bmatrix} \]

应用：邻接矩阵的幂在图上的意义

我们记一个无边权无向图的邻接矩阵为 \(A\)，\(A_{i,j}=1\) 表示两点间有一条边相连，反之则没有。

我们在Floyd算法中有这么一个松弛操作：

for (int k = 1; k <= n; ++k)
    for (int i = 1; i <= n; ++i)
        for (int j = 1; j <= n; ++j)
            if (d[i][j] > d[i][k] + d[k][j])
                d[i][j] = d[i][k] + d[k][j];

这是为了寻找中继点，尝试用来更新两点间的路径。我们惊奇发现，这个流程似乎和矩阵乘法十分相似。

如果我们将他稍微修改一下，变成这样：

int d[N][N], D[N][N];
for (int i = 1; i <= n; ++i)
    for (int j = 1; j <= n; ++j)
        for (int k = 1; k <= n; ++k)
            D[i][j] += d[i][k] * d[k][j];

我们稍加猜想，便可想到 \(D_{i,j}\) 存储的是经过两步，从点i走到点j的方案数。有趣的是，\(D=d^2\) 。实际上，\(d_{i,j}\) 也可以理解为经过一步，从点i走到点j的方案数。

而我们计算 \(D^2=d^4\)，就可以得到经过四步的方案数，\(D*d=d^3\) 则是经过三步的方案数，诸如此类（例如\(d^x*d^y\)，我们枚举 \(k\)，那么 \(d^x_{i,k}*d^y_{k,j}\) 便是走 \(x\) 步从 i 到 k，再走 \(y\) 步从 k到 j的方案数，累加上去后，便得到了走 \(x+y\) 步以从 i 到 j的方案数）。

实际上，结合Floyd算法和矩阵乘法的思想，我们不难得出结论：记 \(A\) 为邻接矩阵，则 \(A^k\) 则记录了某个点走 \(k\) 步后到达另一个点的总方案个数。

有一个机器人在一张无边权无向图（n 点 m 边）的起点1上。每过一秒钟，他便可以进行下面三个操作之一：

停在原地不动

走向一个相邻的点

自爆

问经过时间 \(t\) 后，机器人的行为方案总数有多少？

\(1\leq n \leq 30,0\leq m \leq 100,1\leq t\leq 10^6\)

题目来源：P3758 [TJOI2017]可乐

如果只有第二个操作，那么就是一个板子题了：直接计算邻接矩阵的 t 次方，记为 \(D\)，然后统计 \(\sum\limits_{i=1}^{n}D_{1,i}\) 即可得出方案总数。

对于第一个操作，我们给每个点连上一个自环即可。

对于第三个操作，我们可以建一个新点，来表示这个死亡状态：每个点都向它连边，代表机器人在任意位置都可以自爆；这个点不向外连边，代表机器人自爆之后不可能复活，前往别的位置啥的。（但是自己要往自己连一条边，不然的话所花时间就没法保证完全为 \(t\) 了，不懂的话可以拿样例试一试）

#include<bits/stdc++.h>
using namespace std;
#define LL long long
const int N = 40;
const LL mod = 2017;
struct Matrix {
    int n;
    LL a[N][N];
    void init(int _n, int opt) {
        memset(a, 0, sizeof(a));
            n = _n;
            if (opt == 1)
                for (int i = 1; i <= n; ++i)
                    a[i][i] = 1;
    }
    Matrix operator = (const Matrix &B) {
        this->n = B.n;
        for (int i = 1; i <= n; ++i)
            for (int j = 1; j <= n; ++j)
                this->a[i][j] = B.a[i][j];
        return *this;
    }
    Matrix operator * (const Matrix &B) const {
        Matrix res;
        res.init(n, 0);
        for (int i = 1; i <= n; ++i)
            for (int j = 1; j <= n; ++j)
                for (int k = 1; k <= n; ++k) {
                    res.a[i][j] += (this->a[i][k]) * B.a[k][j];
                    res.a[i][j] %= mod;
                }
        return res;
    }
    void output() {
        for (int i = 1; i <= n; ++i) {
            for (int j = 1; j < n; ++j)
                printf("%lld ", a[i][j]);
            printf("%lld\n", a[i][n]);
        }
    }
};
Matrix quickpow(Matrix P, LL n) {
    if (n == 1) return P;
    Matrix res;
    res.init(P.n, 1);
    while (n) {
        if (n & 1) res = res * P;
        n >>= 1;
        P = P * P;
    }
    return res;
}
int main()
{
    Matrix A;
    int n, m, t;
    scanf("%d%d", &n, &m);
    A.init(n + 1, 1);
    for (int i = 1; i <= n; ++i)
        A.a[i][n + 1] = 1;
    for (int i = 1; i <= m; ++i) {
        int u, v;
        scanf("%d%d", &u, &v);
        A.a[u][v] = A.a[v][u] = 1;
    }
    scanf("%d", &t);
    Matrix D = quickpow(A, t);
    LL ans = 0;
    for (int i = 1; i <= n + 1; ++i)
        ans = (ans + D.a[1][i]) % mod;
    printf("%lld", ans);
    return 0;
}

应用：对DP进行加速

动态规划的求解形式有两种：一是记忆化搜索，二是递推。如果一个基于递推的 DP 方程中仅包含加减等基本运算，便可以使用矩阵乘法对其进行优化。

给定一个长度为 \(n\) 的数组，你现在需要将其分块。

小明要求每个块的长度必须在集合 \(S\) 种，小红要求每个块的长度必须在集合 \(T\) 中，问我们有多少种不同的分块方法？

记最大块长为 \(x\)，有 \(1\leq n\leq 10^{18},1\leq |S|,|T|,x\leq 100\)

P5343 【XR-1】分块

我们直接对两个集合求一下交集集合，方法随意（反正范围就这么点大）。

我们不妨记交集为 \(D\)，记 \(f_i\) 表示前 \(i\) 长度的数列的分块种类，那么有 DP 方程

\[f_0=1,f_i=\sum\limits_{v\in D,v\leq i}f_{i-v} \]

如果线性递推，复杂度显然是 \(O(nx)\) 的，直接 T 飞。不过，我们可以用一手矩阵快速幂来对他进行优化。

已知当 \(i\geq x\) 的时候，一整个求和式子都可以算满，所以我们先预处理出 0 到 x-1 的初始向量 \(\begin{bmatrix}f_{x-1}&f_{x-2} & \cdots&f_1&f_0\end{bmatrix}^T\) ，随后我们直接构造递推矩阵即可。其实并不是很难，照着集合 \(D\) 的情况递推即可。我们假设 \(D=\{1,2,5,7\}\)，那么便有 DP 方程：

\[\begin{bmatrix}f_{i}\\f_{i-1} \\ \cdots\\f_{i-5}\\f_{i-6}\end{bmatrix}= \begin{bmatrix}1&1&0&0&1&0&1\\ 1&0&0&0&0&0&0 \\ \vdots&&&\cdots &&&\vdots\\ 0&0&0&0&1&0&0\\ 0&0&0&0&0&1&0 \end{bmatrix} \begin{bmatrix}f_{i-1}\\f_{i-2} \\ \cdots\\f_{i-6}\\f_{i-7}\end{bmatrix} \]

矩阵快速幂进行优化，可以将复杂度降到 \(O(x^3\log n)\)。

#include<bits/stdc++.h>
using namespace std;
#define LL long long
const int N = 110;
const LL mod = 1e9 + 7;
struct Matrix {
    int n;
    LL a[N][N];
    void init(int _n, int opt) {
        memset(a, 0, sizeof(a));
            n = _n;
            if (opt == 1)
                for (int i = 1; i <= n; ++i)
                    a[i][i] = 1;
    }
    Matrix operator = (const Matrix &B) {
        this->n = B.n;
        for (int i = 1; i <= n; ++i)
            for (int j = 1; j <= n; ++j)
                this->a[i][j] = B.a[i][j];
        return *this;
    }
    Matrix operator * (const Matrix &B) const {
        Matrix res;
        res.init(n, 0);
        for (int i = 1; i <= n; ++i)
            for (int j = 1; j <= n; ++j)
                for (int k = 1; k <= n; ++k) {
                    res.a[i][j] += (this->a[i][k]) * B.a[k][j];
                    res.a[i][j] %= mod;
                }
        return res;
    }
};
Matrix quickpow(Matrix P, LL n) {
    if (n == 1) return P;
    Matrix res;
    res.init(P.n, 1);
    while (n) {
        if (n & 1) res = res * P;
        n >>= 1;
        P = P * P;
    }
    return res;
}
namespace getSelf {
    int n, m, vis1[N], vis2[N];
    //范围比较小，可以通过开一个桶来求交集
    void getArr(int &cnt, int *arr) {
        memset(vis1, 0, sizeof(vis1));
        memset(vis2, 0, sizeof(vis2));
        scanf("%d", &n);
        for (int i = 0, x; i < n; ++i)
            scanf("%d", &x), vis1[x] = 1;
        scanf("%d", &m);
        for (int i = 0, x; i < m; ++i)
            scanf("%d", &x), vis2[x] = 1;
        for (int i = 1; i < N; ++i)
            if (vis1[i] && vis2[i]) arr[++cnt] = i;
    }
}
int main()
{
    //读入数据
    LL n;
    cin >> n;
    int m = 0, v[N];
    getSelf::getArr(m, v);
    //预处理DP
    int x = v[m];
    LL f[N]; memset(f, 0, sizeof(f)); f[0] = 1;
    for (int i = 1; i < x; ++i)
        for (int j = 1; j <= m; ++j)
            if (i >= v[j]) f[i] = (f[i] + f[i - v[j]]) % mod;
    //构造转移矩阵
    Matrix A;
    A.init(x, 0);
    for (int i = 1; i <= m; ++i) A.a[1][v[i]] = 1;
    for (int i = 2; i <= x; ++i) A.a[i][i - 1] = 1;
    Matrix B = quickpow(A, n - x + 1);
    //求和输出
    LL ans = 0;
    for (int i = 1; i <= x; ++i)
        ans = (ans + B.a[1][i] * f[x - i]) % mod;
    printf("%lld", ans);
    return 0;
}

posted @ 2021-08-09 22:22 cyhforlight 阅读(68) 评论(0) 编辑收藏举报

会员力量，点亮园子希望

刷新页面返回顶部

cyhforlight

2021 暑假训练学习 矩阵快速幂