浅谈矩阵乘法

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\)为区间求和，线段树都可以完成。

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