模板

在使用的时候一定记得要初始化！！！

#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;
}

应用：优化线性递推式

求斐波那契数列的第 项 。（#card=math&code=f1%3Df_2%3D1%2Cf_n%3Df%7Bn-1%7D%2Bf_%7Bn-2%7D%28n%20%5Cgeq%203%29&id=YH9Lb)）

递推的复杂度是 #card=math&code=O%28n%29&id=TCYlu)，现在不能接受。

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

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

对于 ，有

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

#include<bits/stdc++.h>
using namespace std;
#define LL long long
const int N = 5;
const LL mod = 1e9 + 7;
//略去Matrix和quickpow的部分，以下同理
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;
}

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

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

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

我们记一个无边权无向图的邻接矩阵为 ， 表示两点间有一条边相连，反之则没有。
我们在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];

我们稍加猜想，便可想到 存储的是经过两步，从点i走到点j的方案数。有趣的是， 。实际上， 也可以理解为经过一步，从点i走到点j的方案数。

而我们计算 ，就可以得到经过四步的方案数， 则是经过三步的方案数，诸如此类（例如，我们枚举 ，那么 便是走 步从 i 到 k，再走 步从 k到 j的方案数，累加上去后，便得到了走 步以从 i 到 j的方案数 ）。

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

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

1. 停在原地不动
2. 走向一个相邻的点
3. 自爆

问经过时间 后，机器人的行为方案总数有多少？ 题目来源：P3758 [TJOI2017]可乐

如果只有第二个操作，那么就是一个板子题了：直接计算邻接矩阵的 t 次方，记为 ，然后统计 即可得出方案总数。

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

#include<bits/stdc++.h>
using namespace std;
#define LL long long
const int N = 40;
const LL mod = 2017;
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 方程中仅包含加减等基本运算，便可以使用矩阵乘法对其进行优化。

给定一个长度为 的数组，你现在需要将其分块。 小明要求每个块的长度必须在集合 种，小红要求每个块的长度必须在集合 中，问我们有多少种不同的分块方法？ 记最大块长为 ，有

我们直接对两个集合求一下交集集合，方法随意（反正范围就这么点大）。
我们不妨记交集为 ，记 表示前 长度的数列的分块种类，那么有 DP 方程

如果线性递推，复杂度显然是 #card=math&code=O%28nx%29&id=bjmYf) 的，直接 T 飞。不过，我们可以用一手矩阵快速幂来对他进行优化。
已知当 的时候，一整个求和式子都可以算满，所以我们先预处理出 0 到 x-1 的初始向量 ，随后我们直接构造递推矩阵即可。其实并不是很难，照着集合 的情况递推即可。我们假设 ，那么便有 DP 方程：

矩阵快速幂进行优化，可以将复杂度降到 #card=math&code=O%28x%5E3%5Clog%20n%29&id=SYIHO)。

#include<bits/stdc++.h>
using namespace std;
#define LL long long
const int N = 110;
const LL mod = 1e9 + 7;
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;
}