矩阵快速幂学习笔记

Wed Dec 18 2024

9 分钟

1415 字

这矩阵快速幂多是一件美事啊
只用推式子，而且这式子也是十分好推。

P1962 斐波那契数列 #

是个人第一眼看到这题，第一眼都会想着递推。但是如果真这么简单的话，它就不会是绿题了。我们一看数据范围： $n \le 2^{63}$ 。 $\mathcal{O}(n)$ 不给你全部炸完？
先来推式子： $f_{i} = f_{i - 1} + f_{i - 2} , f_{i - 1} = f_{i - 1}$ 。我们把他整理一下：

\left\{ \begin{aligned} f_{i} & = 1 \times f_{i - 1} + 1 \times f_{i - 2} \\ f_{i - 1} & = 1 \times f_{i - 1} + 0 \times f_{i - 2} \end{aligned} \right.

我们一看，这可以用矩阵来表示啊：

\begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix} \begin{bmatrix} f_{i - 1} \\ f_{i - 2} \end{bmatrix} = \begin{bmatrix} f_{i} \\ f_{i - 1} \end{bmatrix}

而我们又看：

\begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix}^1 \begin{bmatrix} f_{i - 1} \\ f_{i - 2} \end{bmatrix} = \begin{bmatrix} f_{i} \\ f_{i - 1} \end{bmatrix} \\ \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix}^2 \begin{bmatrix} f_{i - 2} \\ f_{i - 3} \end{bmatrix} = \begin{bmatrix} f_{i} \\ f_{i - 1} \end{bmatrix}

发现什么端倪没有？如果我们一直往前推，定能推到 $\begin{bmatrix}f_{2} \\ f_{1}\end{bmatrix}$ 。而我们看：往前推一次，要乘转矩阵的 $1$ 次方， $f_{i}$ 变成 $f_{i - 1}$ 。而我们最终要从 $f_{n}$ 变成 $f_{2}$ ，也就是 $f_{n - (n - 2)}$ ，那么我们就可以知道，推到 $\begin{bmatrix}f_{2} \\ f_{1}\end{bmatrix}$ 要用原矩阵 $\begin{bmatrix}1 \\ 1 \end{bmatrix}$ 乘转移矩阵 $\begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix}$ 的 $n - 2$ 次方。也就是 $\begin{bmatrix}f_{n} \\ f_{n - 1}\end{bmatrix} = \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix}^{n - 2}\begin{bmatrix}1 \\ 1 \end{bmatrix}$ 。那么我们就可以用矩阵快速幂来做咯。注意：矩阵乘法没有交换律，一定注意矩阵乘法的顺序。

code:

CPP

#include<bits/extc++.h>
#define int long long
using namespace std;
const int mod = 1e9 + 7;
int n;
struct mat{
    int a[2][2];
    mat(int x = 0){memset(a,0,sizeof a);a[0][0] = a[1][1] = x;};
    mat operator*(mat x)
    {
        mat ret;
        for (int i = 0; i < 2; i++)
            for (int j = 0; j < 2; j++)
                for (int k = 0; k < 2; k++)
                    ret.a[i][j] = (ret.a[i][j] + a[i][k] * x.a[k][j] % mod) % mod;
        return ret;
    }
}base,tmp;
template<typename type>
type binpow(type x,int y)
{
    type ret(1);
    while (y)
    {
        if (y & 1)
            ret = ret * x;
        x = x * x;
        y >>= 1;
    }
    return ret;
}
signed main()
{
    scanf("%lld",&n);
    if (n == 1)
    {
        printf("1");
        return 0;
    }
    base.a[0][0] = 1;
    base.a[0][1] = 1;
    base.a[1][0] = 1;
    base.a[1][1] = 0;
    base = binpow(base,n - 2);
    tmp.a[0][0] = tmp.a[0][1] = 1;
    printf("%lld",(tmp * base).a[0][0]);
    return 0;
}

P6569 [NOI Online #3 提高组] 魔法值 #

我们看： $n \le 100$ ，用邻接矩阵存图吧，在邻接矩阵中， $i$ 到 $j$ 有边， $dis_{i,j}$ 就为 $1$ ，反之为 $0$ 。然后我们就可以写出这样的一个柿子：

f_{i,j} = \bigoplus\limits_{k = 1}^{n}f_{j - 1,k} \times dis_{k,i}

我们看，这玩意，似乎有点像什么？矩阵乘法啊。只不过是把求和换成了异或而已。然后，我们就可以定义新的矩阵乘法了。但是， $f_{i,j}$ 和式子里的顺序似乎反了，于是我们直接把 $i,j$ 的定义反一下：

f_{i,j} = \bigoplus\limits_{k = 1}^{n}f_{i - 1,k} \times dis_{k,j}

也是十分顺眼。

然后，我们看，矩阵乘法，那不得用快速幂啊，但是快速幂要用到乘法结合律，我们来看看这玩意有没有结合律。啊，好像没有。但是不用快速幂又没法做啊，硬着头皮证一下吧。
发现 $dis$ 一定是一个 01 矩阵，似乎这个特殊性质有点用？
我们要证结合律，那么就是要证 $(A \times B) \times C = A \times (B \times C)$ 。根据定义， $(A \times B) \times C$ 的 $x,y$ 处的元素是 $\bigoplus\limits_{i = 1}^{n}(\bigoplus\limits_{j = 1}^{n}A_{x,j} \times B_{j,i}) \times C_{i,y}$ 。但是异或是没有分配律的，所以这括号现在还不能拆，得三思而后行。我们看 $C_{i,y}$ 只有 $0$ 和 $1$ 两种情况，分类讨论一下？

当 $C_{i,y} = 0$ 时，那一大坨都给乘没了，所以可以拆。
当 $C_{i,y} = 1$ 时，乘完还是那坨，也可以。

综上所述，那个括号是可以拆的，当且仅当 $C$ 是 01 矩阵时。巧了，我们的邻接矩阵就是 01 矩阵，这不爽了？

根据上面推的式子： $f_{i} = f_{i - 1} \times dis$ ，其中 $f_{i}$ 与 $f_{i - 1}$ 都是 $1$ 行 $n$ 列的矩阵。那么显然有： $f_{i} = f_{0} \times dis^{i}$ 。这不就是快速幂吗？但是，我们看，查询有 $100$ 个，而矩阵乘法是 $\mathcal{O}(n^3)$ 的，加上快速幂的总复杂度就是 $\mathcal{O}(qn^3\log a_i)$ ，似乎不大行，于是我们考虑预处理。把 $dis$ 的 $2^{x}$ 都预处理出来，查询的时候循坏利用，就可以了。

code:

CPP

#include<bits/extc++.h>
#define int long long
using namespace std;
int n,m,q;
int f[105];
struct mat
{
    int r,c;
    int a[105][105];
    mat(){memset(a,0,sizeof a);};
    mat(int x,int r,int c):r(r),c(c)
    {
        memset(a,0,sizeof a);
        for (int i = 1; i <= 100; i++)
            a[i][i] = x;
    };
    int *operator[](int x){return a[x];}
    friend mat operator*(mat x,mat y)
    {
        mat ret(0,x.r,y.c);
        for (int i = 1; i <= x.r; i++)
            for (int j = 1; j <= y.c; j++)
                for (int k = 1; k <= x.c; k++)
                    ret[i][j] ^= x[i][k] * y[k][j];//异或矩阵乘法
        return ret;
    }
}dis[64];
signed main()
{
    cin >> n >> m >> q;
    for (int i = 1; i <= n; i++)
        cin >> f[i];
    dis[0].r = dis[0].c = n;
    int u,v;
    for (int i = 1; i <= m; i++)
    {
        cin >> u >> v;
        dis[0][u][v] = dis[0][v][u] = 1;
    }
    for (int i = 1; i < 64; i++)
        dis[i] = dis[i - 1] * dis[i - 1];//预处理dis的2^x次方
    int x;
    for (int i = 1; i <= q; i++)
    {
        cin >> x;
        mat origin(0,n,1);
        for (int j = 1; j <= n; j++)
            origin[j][1] = f[j];//f[0]初始化，也就是初始化原矩阵
            //本人习惯打的n行1列的矩阵，但是上面的式子推着推着就成了1行n列，所以就有区别
        for (int j = 0; j < 64; j++)
            if ((x >> j) & 1)
                origin = dis[j] * origin;
        cout << origin[1][1] << endl;
    }
    return 0;
}