浅谈矩阵乘法

1.矩阵乘法

给出两个矩阵$A,B$，令$C$为$A\times B$的结果，若$A$是一个$n\times m$的矩阵，$B$是一个$m\times p$的矩阵，则$C$是一个$n\times p$的矩阵。其中有：

$$\forall i\in [1,n],j\in [1,p],C_{i,j}=\sum\limits_{k=1}^{m} A_{i,k}\times B_{k,j}$$

注意到由于矩阵乘法的定义，只有当$A$的列数等于$B$的行数时才可以相乘，否则无意义。

柿子的定义可能比较抽象，难以理解，那我们就举个例子：

$$\begin{bmatrix}1&2&3\\4&5&6\end{bmatrix} \cdot\begin{bmatrix}1&2\\3&4\\5&6\end{bmatrix}$$

根据柿子手摸可以得到结果：

$$\begin{bmatrix}22&28\\49&64\end{bmatrix}$$

形象化地说，$C_{i,j}$就是把$A$的第$i$行和$B$的第$j$列的$m$个数一个个乘起来再相加。

参考代码：

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

class Matrix {
   public:
    int n, m;
    int val[55][55];
    Matrix() { memset(val, 0, sizeof(val)); }
    void init() {
        for (int i = 1; i <= 50; ++i) {
            for (int j = 1; j <= 50; ++j) {
                val[i][j] = 1;
            }
        }
    }
};

Matrix operator*(const Matrix &a, const Matrix &b) {
    int i, j, k;
    Matrix res;
    for (k = 1; k <= a.m; ++k) {
        for (i = 1; i <= a.n; ++i) {
            for (j = 1; j <= b.m; ++j) {
                res.val[i][j] += a.val[i][k] * b.val[k][j];
            }
        }
    }

    return res;
}

signed main() {
    Matrix a, b, ans;
    int n, m, p;
    int i, j;

    scanf("%d%d%d", &n, &m, &p);
    a.n = n, a.m = m;
    b.n = m, b.m = p;
    for (i = 1; i <= n; ++i) {
        for (j = 1; j <= m; ++j) {
            scanf("%d", &a.val[i][j]);
        }
    }

    for (i = 1; i <= m; ++i) {
        for (j = 1; j <= p; ++j) {
            scanf("%d", &b.val[i][j]);
        }
    }

    ans = a * b;
    for (i = 1; i <= n; ++i) {
        for (j = 1; j <= p; ++j) {
            printf("%d\n", ans.val[i][j]);
        }
    }
    return 0;
}

2.矩阵快速幂

矩阵乘法满足结合律，即对于矩阵$A,B,C$，有

$$(A\times B)\times C = A \times (B\times C)$$

所以我们可以对于矩阵进行经典快速幂运算。

这里要注意，矩阵乘法不满足交换律，即$A\times B \not=B\times A$

洛谷P3390【模板】矩阵快速幂

参考代码：

#include <bits/stdc++.h>
using namespace std;
const int p = 1e9 + 7;
#define int long long

inline int read()
{
    int ans = 0, f = 1;
    char ch = getchar();
    while (ch > '9' || ch < '0') {
        if (ch == '-')
            f = -1;
        ch = getchar();
    }
    while (ch >= '0' && ch <= '9') {
        ans = (ans << 1) + (ans << 3) + (ch ^ '0');
        ch = getchar();
    }
    return ans * f;
}

int n, k;

struct Matrix {
    int a[105][105];
    Matrix()
    {
        memset(a, 0, sizeof a);
    }
    inline void init()
    {
        for (int i = 1; i <= n; ++i)
            a[i][i] = 1;
    }
} t, ans;

Matrix operator*(const Matrix& x, const Matrix& y)
{
    Matrix z;
    for (int k = 1; k <= n; ++k)
        for (int i = 1; i <= n; ++i)
            for (int j = 1; j <= n; ++j)
                z.a[i][j] = (z.a[i][j] + x.a[i][k] * y.a[k][j] % p) % p;
    return z;
}

Matrix power(Matrix a, int b)
{
    Matrix ans, base = a;
    ans.init();
    while (b) {
        if (b & 1)
            ans = ans * base;
        base = base * base;
        b >>= 1;
    }
    return ans;
}

signed main()
{
    n = read(), k = read();
    register int i, j;
    ans.init();

    for (i = 1; i <= n; ++i) {
        for (j = 1; j <= n; ++j) {
            t.a[i][j] = read();
        }
    }
    ans = power(t, k);

    for (i = 1; i <= n; ++i) {
        for (j = 1; j <= n; ++j)
            printf("%lld ", ans.a[i][j]);
        puts("");
    }

    return 0;
}

3.应用

我们来看一个经典例题：

给出$n$，求斐波那契数列的第$n$项mod 998244353
$1\le n \le 10^9$

暴力递推显然满足不了$10^9$的数据范围，所以我们需要加速这个过程。

加速的方法很多，比如特征方程求通项公式等数学方法；这里给出矩阵快速幂的解法。

我们构造一个矩阵：

$$\begin{bmatrix} f_i & f_{i-1} \end{bmatrix}$$

其中$f_i$表示斐波那契数列的第$i$项。

我们把这个矩阵乘上一个base矩阵：

$$\begin{bmatrix}1 & 1\\1& 0\end{bmatrix}$$

得到：

$$\begin{bmatrix}f_i\times 1+f_{i-1}\times 1&f_i\times 1 + f_{i-1}\times 0\end{bmatrix}$$

即：

$$\begin{bmatrix}f_{i+1}&f_i\end{bmatrix}$$

发现什么？我们通过一次矩阵乘法的运算即得到了斐波那契数列的下一项，完成了一次递推。同时，矩阵可以进行快速幂，所以我们可以通过乘上base矩阵的$n$次幂来快速递推，复杂度被优化到了$O(\log n)$

到这里可以发现，矩阵快速幂可以加速递推式的计算过程。

我们再次观察上述求解过程：到底是怎么构造出原矩阵和base矩阵的呢？

斐波那契数列的通项公式为：$f_i=f_{i-1}+f_{i-2}$，所以我们要知道$f_{i+1}$必须要知道$f_i$和$f_{i-1}$。因为矩阵乘法是两个矩阵内部的运算，不可能凭空出现一个新的值，所以我们肯定要在原矩阵里面记录下这两个值才能得到$f_{i+1}$。

所以现在我们有矩阵$\begin{bmatrix}f_i & f_{i-1}\end{bmatrix}$，要乘上一个base矩阵之后得到$\begin{bmatrix}f_{i+1} & f_{i}\end{bmatrix}$，又

$$f_{i+1}= f_i \times 1+f_{i-1}\times 1\\f_i=f_i\times 1 + f_{i-1}\times 0$$

将系数联立起来即可得到base矩阵。

4.例题

1.洛谷P1939【模板】矩阵加速（数列）

给出$n$，已知$f_i=f_{i-1}+f_{i-3}$，求$f_n$
$1\le n\le 10^9$

我们构造矩阵$\begin{bmatrix}f_i&f_{i-1}&f_{i-2}\end{bmatrix}$，希望得到结果为$\begin{bmatrix}f_{i+1}&f_i&f_{i-1}\end{bmatrix}$

显然有：

$$f_{i+1}=f_i\times 1+f_{i-1}\times 0+f_{i-2}\times 1\\ f_i=f_i\times 1 + f_{i-1}\times 0+f_{i-2}\times 0\\ f_{i-1}=f_i\times 0 + f_{i-1}\times 1 +f_{i-2}\times 0$$

所以构造出base矩阵：

$$\begin{bmatrix}1&1&0\\0&0&1\\1&0&0\end{bmatrix}$$

2.[NOI2012] 随机数生成器

给出$n,a,c$，已知$f_i=af_{i-1}+c$，求$f_n$
$1\le n \le 10^{18}$

这里发现$f_i$只与$f_{i-1}$有关，所以我们构造出原始矩阵：

$$\begin{bmatrix}f_i& c\end{bmatrix}$$

显然有：

$$f_{i+1}=a\times f_i+1\times c\\ c=0\times f_i+1\times c$$

所以我们乘上：

$$\begin{bmatrix}a&0\\1&1\end{bmatrix}$$

即可得到：

$$\begin{bmatrix}f_{i+1}&c\end{bmatrix}$$

这道题目原题好像裸的矩阵快速幂会炸longlong要加上龟速乘，不过这跟无辜的矩阵有什么关系呢

3.LOJ#10222 佳佳的 Fibonacci

设$f_i$为斐波那契数列的第$i$项，$g_i= \sum\limits_{j=1}^{i}f_j\times j$，求$g_n$
$1\le n\le 2^{31}-1$

我们发现

$$g_i=g_{i-1}+i\times f_i\\ f_i=f_{i-1}+f_{i-2}$$

所以我们可以根据这个递推式构造矩阵：

$$\begin{bmatrix}g_i&if_i&(i-1)f_{i-1}&f_i&f_{i-1}\end{bmatrix}$$

因为

$$g_{i+1}=g_i+(i+1)f_{i+1}\\ (i+1)f_{i+1}=if_i+f_i+(i-1)f_{i-1}+2f_{i-1}\\ if_i=if_i\\ f_{i+1}=f_i+f_{i-1}\\ f_i=f_i$$

所以我们得到base矩阵为

$$\begin{bmatrix}1&0&0&0&0\\1&1&1&0&0\\1&1&0&0&0\\1&1&0&1&1\\2&2&0&1&0\end{bmatrix}$$

将原矩阵乘上base矩阵即可得到：

$$\begin{bmatrix}g_{i+1}&(i+1)f_{i+1}&if_i&f_{i+1}&f_i\end{bmatrix}$$

4.[HNOI2011]数学作业

定义$f_i$为从$1$到$i$的所有整数写在一起得到的数，如$f_7=1234567,f_{13}=12345678910111213$，求$f_n\mod m$
$1\le n\le 10^{18}$

显然有$f_i=f_{i-1}\times 10^{len(i)}+i$，其中$len(i)$表示$i$的位数，所以可以就此进行矩阵加速递推。

原矩阵：

$$\begin{bmatrix}f_i&i&1\end{bmatrix}$$

base矩阵：

$$\begin{bmatrix}10^{len(i+1)}&0&0\\1&1&0\\1&1&1\end{bmatrix}$$

但是我们发现$len(i)$不好维护，考虑到$1\le n\le 10^{18}$，所以$len$的值只有18种左右，枚举位数分开矩乘即可。

5.[SNOI2017]礼物

设$f_1=1,s_i=\sum\limits_{j=1}^{i}f_j,f_i=s_{i-1}+i^k$，求$f_n$
$1\le n\le 10^{18},1\le k\le 10$

先推推柿子：

$$s_i=s_{i-1}+f_i=2s_{i-1}+i^k$$

这样我们得到了$s$的递推式，尝试矩阵优化；但是我们发现不好直接维护$i^k$，考虑二项式定理：

$$(x+1)^k=\sum\limits_{i=0}^{k}\begin{pmatrix}k\\i\end{pmatrix}x^i$$

于是便可以构造矩阵：

$$\begin{bmatrix}s_i&i^k&i^{k-1}&\cdots& i^0\end{bmatrix}$$

我们希望得到

$$\begin{bmatrix}s_{i+1}&(i+1)^k&(i+1)^{k-1}&\cdots& (i+1)^0\end{bmatrix}$$

于是便可以推出base矩阵：

$$\begin{bmatrix}2&0&0&0&0&\cdots &0\\\begin{pmatrix}k\\0\end{pmatrix}&\begin{pmatrix}k\\0\end{pmatrix}&0&0&0&\cdots &0\\ \begin{pmatrix}k\\1\end{pmatrix}&\begin{pmatrix}k\\1\end{pmatrix}&\begin{pmatrix}k-1\\0\end{pmatrix}&0&0&\cdots &0\\ \begin{pmatrix}k\\2\end{pmatrix}&\begin{pmatrix}k\\2\end{pmatrix}&\begin{pmatrix}k-1\\1\end{pmatrix}&\begin{pmatrix}k-2\\0\end{pmatrix}&0&\cdots &0\\ \vdots&\vdots&\vdots&\vdots&\vdots&\ddots&\vdots\\ \begin{pmatrix}k\\k\end{pmatrix}&\begin{pmatrix}k\\k\end{pmatrix}&\begin{pmatrix}k-1\\k-1\end{pmatrix}&\begin{pmatrix}k-2\\k-2\end{pmatrix}&\begin{pmatrix}k-3\\k-3\end{pmatrix}&\cdots&\begin{pmatrix}0\\0\end{pmatrix} \end{bmatrix}$$

这里定义$\begin{pmatrix}0\\0\end{pmatrix}=1$

注意到$k$非常小，所以以上做法可以通过。

6.[THUSCH2017] 大魔法师

给出$n$个元素，每个元素有$3$个值：$A,B,C$，共进行$m$次操作，每次选取一个区间$[l,r]$进行如下几种操作之一：
$1.A_i\leftarrow A_i+B_i$
$2.B_i\leftarrow B_i+C_i$
$3.C_i\leftarrow C_i+A_i$
$4.A_i\leftarrow A_i+v$
$5.B_i\leftarrow B_i\times v$
$6.C_i\leftarrow v$
$7.$求$\sum\limits_{i=l}^{r}A_i,\sum\limits_{i=l}^{r}B_i,\sum\limits_{i=l}^{r}C_i$
请对于每个操作$7$输出结果
$1\le n,m\le 2.5\times 10^5$

我们把每个元素看成一个矩阵：

$$\begin{bmatrix}A_i&B_i&C_i&1\end{bmatrix}$$

显然每次操作可以用矩阵乘法来表示，例如操作$1$的base为：

$$\begin{bmatrix}1&0&0&0\\1&1&0&0\\0&0&1&0\\0&0&0&1\end{bmatrix}$$

操作$2-6$同理

我们对于这$n$个矩阵开一棵线段树，显然操作$1-6$即为区间乘，操作$7$为区间求和，线段树都可以完成。

注意到矩阵存在分配律，所以可以用懒标记实现。