算法矩阵快速幂

方今之时，臣以神遇而不以目视，官知止而神欲行，依乎天理，批大郤，导大窾，因其固然，技经肯綮之未尝，而况大軱乎！

《庄子·养生主》

昔庖丁历十九岁，所解数千牛后方有依乎天理之技。解牛固然，代码亦固然。千题者与十题者孰轻重，其安有疑？所以分数等第之高低者，实乃其能力之高低上下也。所以致其能力之高低上下者，无乃手熟与手生乎？

矩阵快速幂，是一种可以对于特殊的动态规划进行时间优化的算法。

理论推导

向量、矩阵

所谓向量，一般可以写成一个由若干个数组成的一个序列，此若干个数可以横写，也可以竖写，其表示的意义相同，如右侧和下侧。在算法上，一般用向量表示一个动态规划中的状态。

\[A=\begin{bmatrix} a_1 &a_2 &a_3 &\cdots &a_n\end{bmatrix}\]

\[A=\begin{bmatrix} a_1 \\ a_2 \\ a_3 \\ \cdots \\ a_n\end{bmatrix}\]

矩阵是一个数按行列排列而成的一种数据单位，一般一个 \(n\) 行 \(m\) 列的矩阵我们成为 \(n\times m\) 的矩阵。特别地，对于 \(n=1\) 或 \(m=1\) 的情况我们也称之为 \(n\times 1\) 或 \(1\times m\) 的向量。在矩阵中，第 \(i\) 行第 \(j\) 列的元素下标为 \(ij\) ，如下图。在算法中，矩阵一般表示两个向量表示的状态之间的转化规则，即动态规划的转移方程。

\[S=\begin{bmatrix} s_{11} &s_{12} &s_{13} &\cdots &s_{1m} \\ s_{21} &s_{22} &s_{23} &\cdots &s_{2m} \\ s_{31} &s_{32} &s_{33} &\cdots &s_{3m} \\ \vdots &\vdots &\vdots &\ddots &\vdots \\ s_{n1} &s_{n2} &s_{n3} &\cdots &s_{nm}\end{bmatrix}\]

矩阵乘法

从算法角度讲，矩阵乘法就是对于动态规划转移状态的计算。关于矩阵乘法的解释，在此引入一个例子辅助理解。

假若我们已知向量 \(A_1=\begin{bmatrix} a_1 &a_2\end{bmatrix}\) ，想要通过一个矩阵运算得到 \(A_2=\begin{bmatrix} a_1+2a_2 &3a_1-4a_2\end{bmatrix}\) 。在这里向量 \(A_2\) 的第一项等于向量 \(A_1\) 的第一项加上二倍第二项，向量 \(A_2\) 的第二项等于向量 \(A_1\) 的三倍第一项减去四倍第二项，故我们可以将其总结为一个转移矩阵 \(S\)，并列出等式：

\[A_2=SA_1\]

\[\begin{bmatrix} a_1+2a_2 \\ 3 a_1-4a_2\end{bmatrix}=\begin{bmatrix} 1 &2 \\ 3 &-4 \\\end{bmatrix}\begin{bmatrix} a_1 \\ a_2\end{bmatrix}\]

矩阵 \(S\) 中的第一行的两项表示运算结果 \(A_2\) 的第一项由 \(A_1\) 的两项按照矩阵中的权值相乘后相加而得，第二行亦然。

同时，从数学上可以总结出一个规律：若是一个 \(n \times m\) 与一个 \(m \times p\) 的矩阵相乘，为了保持系数相匹配，答案便是一个 \(n \times p\) 的矩阵。注意：矩阵乘法不具有交换律，作乘法时须保证前项的列数等于后项的行数。

正因此，我们一般在做矩阵乘法时将向量竖写，并作为乘法的后项。这样根据乘法规则，相乘的结果向量也是竖写的。

通过构建转移矩阵，便可以用矩阵乘法一定程度上代替原来的动态规划。

单位矩阵

与数字乘法类似，矩阵也有单位矩阵。

单位矩阵一定是方阵，主对角线，即从左上到右下的对角线上的数均是 \(1\) ，其他部分为 \(0\) 。如此的单位矩阵可以保证无论作为乘法的前项还是后项，乘法运算之后的结果仍然是原矩阵。

同时，任何方阵的零次方都是对应大小的单位矩阵。

特殊的动态规划

矩阵乘法可以代替一类特殊的动态规划。

考虑这样一道典型的题目：求斐波那契数列的第 \(n\) 项。

考虑传统的动态规划做法，显然若是设 \(F(n)\) 为数列的第 \(n\) 项，则显然有动态规划公式 \(F(n)=F(n-1)+F(n-2)\) ，可以在 \(O(n)\) 的时间内求出第 \(n\) 项。

但是，这个动态规划具有一个特点： \(F(n)\) 的取值是由第 \(F(n-1)\) 和 \(F(n-2)\) 项的线性相加，故若是已知了 \(F(n-1)\) 和 \(F(n-2)\) ，可以推出来 \(F(n)\) 和 \(F(n-1)\) 。可以试图构建这样一个动态规划状态向量：设 \(A_n=\begin{bmatrix} F(n) &F(n-1)\end{bmatrix}\) ，那么在 \(A_{n+1}\) 与 \(A_n\) 之间便可以建立转移矩阵 \(S\) ，使 \(S\) 满足斐波那契数列的规律。如下：

\[A_{n+1}=SA_n\]

\[\begin{bmatrix} F(n+1) \\ F(n)\end{bmatrix}=\begin{bmatrix} F(n)+F(n-1) \\ F(n)\end{bmatrix}=\begin{bmatrix} 1 &1 \\ 1 &0 \\\end{bmatrix}\begin{bmatrix} F(n) \\ F(n-1)\end{bmatrix}\]

这样，便成功地用矩阵乘法来表示了原来的动态规划，可以使用接下来的矩阵快速幂。

此类特殊的动态规划的特点是新状态的各个量都是原状态下量的线性组合，即新状态下的每个量都是由每个原状态量按照固定的权值相加形成。其在矩阵乘法下的意义就是，转移矩阵中的每个量都是常量。在这种状况下，矩阵快速幂才能够实现。

矩阵快速幂

承接上部分的斐波那契数列的例子，注意例中的 \(S\) 是方阵，可以自乘，故由 \(A_n=SA_{n-1}\) 可以推出

\[A_n=SA_{n-1}=S^2A_{n-2}=\cdots=S^{n-2}A_2\]

数学上可以证明，矩阵相乘有结合律，故可以先算 \(S^{n-2}\) ，再与 \(A_2\) 相乘。这个计算 \(S^{n-2}\) 的过程可以利用矩阵快速幂解决，其计算过程与数字快速幂并无太大区别，具体可见代码实现。正因如此，才能够实现以 \(O(\log n)\) 级别实现对时间的优化，胜过一般动态规划。

代码实现

矩阵相关

在此定义方阵结构体，并以重载运算符的方式定义矩阵乘法。

struct matrix
{
    static const int l = 2;
    long long val[l][l];

    matrix();
    matrix operator*(const matrix &) const;
    void unit();
};

其中的静态常量 l 是根据题意定义的方阵边长上限。在运算中，下标均从 \(0\) 开始算起。

对于结构体的构造函数，直接内存清零即可。

matrix::matrix()
{
    memset(val, 0, sizeof(val));
}

对于矩阵的单位化，即转化为单位矩阵，先清零，再把主对角线上的数设为 \(1\) 。

void matrix::unit()
{
    memset(val, 0, sizeof(val));
    for(register int i = 0; i < l; i++)
        val[i][i] = 1;
}

而至于乘法的重载，考虑循环对输出矩阵的第 i 行第 j 列进行计算，其值等于乘法左项矩阵的第 i 行的每个数与乘法右项矩阵第 j 列的每个数分别相乘的结果再相加。语言解释可能略显繁杂，代码显得更加清晰。

matrix matrix::operator*(const matrix &mat) const
{
    matrix res;
    for (register int i = 0; i < l; i++)
        for (register int j = 0; j < l; j++)
            for (register int k = 0; k < l; k++)
                res.val[i][j] += val[i][k] * mat.val[k][j];
    return res;
}

而对于一般的题目，避免数据过大一般会给予模数。设模数为 MOD ，则矩阵乘法可改写为以下代码。

matrix matrix::operator*(const matrix &mat) const
{
    matrix res;
    for (register int i = 0; i < l; i++)
        for (register int j = 0; j < l; j++)
            for (register int k = 0; k < l; k++)
                res.val[i][j] = (res.val[i][j] + val[i][k] * mat.val[k][j]) % MOD;
    return res;
}

由此可知，单次方阵乘法的复杂度是 \(O(l^3)\) ，一个 \(n \times m\) 和 \(m \times p\) 的矩阵相乘，乘法复杂度是 \(O(nmp)\) 。

快速幂代码

函数中的形参 a 是待运算的矩阵， x 是次数。整体和一般的快速幂代码实现几乎相同。

matrix ksm(matrix a, long long x)
{
    matrix ans = a;
    x--;
    for(;x ;x >>= 1)
    {
        if(x & 1)
            ans = ans * a;
        a = a * a;
    }
    return ans;
}

由于考虑到矩阵快速幂和一般快速幂的写法几乎相同，考虑用 C++ 泛型进行模板化，使其同时能够进行矩阵运算和数字运算。

template <typename T>
T ksm(T a, long long x)
{
    T ans = a;
    x--;
    for (; x; x >>= 1)
    {
        if (x & 1)
            ans = ans * a;
        a = a * a;
    }
    return ans;
}

函数返回运算结果。此函数对于 \(x \geqslant 1\) 的情况均能返回正确结果。

模板与例题

P3390 【模板】矩阵快速幂 - 洛谷计算机科学教育新生态
 题解洛谷P3390 【模板】矩阵快速幂 – MindStation

P1939 【模板】矩阵加速（数列） - 洛谷计算机科学教育新生态

P1962 斐波那契数列 - 洛谷计算机科学教育新生态
 题解洛谷P1962 斐波那契数列 – MindStation

P1349 广义斐波那契数列 - 洛谷计算机科学教育新生态

P2044 [NOI2012]随机数生成器 - 洛谷计算机科学教育新生态

P2151 [SDOI2009]HH去散步 - 洛谷计算机科学教育新生态

P6772 [NOI2020]美食家 - 洛谷计算机科学教育新生态

P6569 [NOI Online #3 提高组]魔法值 - 洛谷计算机科学教育新生态

参考资料

俞华程《矩阵乘法在信息学中的应用》

2020-8-3 Update: 增加了模板与例题。
2020-8-4 Update: 优化了快速幂函数的代码，增加了单位矩阵的介绍及代码实现。
2020-8-27 Update: 增加例题（NOI2020 美食家）
2020-9-10 Update: 增加例题（NOI Online#3 魔法值）
2020-9-10 Update: 修改 ksm() 函数，使泛型支持整数运算。