线性代数 - 矩阵快速幂相关

一、加速递推

Ex.1 P1349 广义斐波那契数列

1. 题面梗概

$a_{n} = p \times a_{n-1} + q \times a_{n-2}$，给定 $a_{1},a_{2},p,q,n,m$，求 $a_{n} \bmod m$。

$n \leq 2^{31}-1$

2. 逻辑推导

这就是很经典的矩阵快速幂加速递推。

由于矩阵可以将一个向量进行线性变换，成为另外一个向量。因此可以加速具有明显线性关系的递推式。

因为每次变换只涉及到前两次的状态，所以我们需要一个矩阵将 $\begin{bmatrix}F_{n-1} \\ F_{n-2}\end{bmatrix}$ 变化为 $\begin{bmatrix} F_{n} \\ F_{n-1}\end{bmatrix}$，接着我们采用待定系数法，设转移矩阵为 $\begin{bmatrix}a & b \\ c & d\end{bmatrix}$，那么：

$$\begin{bmatrix} F_{n} \\ F_{n-1} \end{bmatrix} = \begin{bmatrix} a & b \\ c & d \\ \end{bmatrix} \cdot \begin{bmatrix} F_{n-1} \\ F_{n-2} \\ \end{bmatrix}$$

很容易能解出，这个转移矩阵是 $T = \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix}$。将初始向量进行一次变换，可以得到 $\begin{bmatrix}F_{3} \\ F_{2}\end{bmatrix} = T \begin{bmatrix}1 \\ 1 \end{bmatrix}$，接着根据数学归纳法，可以得到：$\begin{bmatrix}F_{n} \\ F_{n-1}\end{bmatrix} = T^{n-2} \begin{bmatrix}1 \\ 1\end{bmatrix}$。

由于矩阵乘法满足结合律，所以可以将矩阵的幂乘使用快速幂的方式来处理。

3. 代码实现

这是一份简易的矩阵快速幂应用，可以视为模版扩展。

▶

AcCode

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

const int maxsz = 100 + 5;

ll n, mod, p, q, a1, a2;

struct Matrix {
    ll m[maxsz][maxsz];
    int sz;
    Matrix( int s = 0 ) : sz( s ) { memset( m, 0, sizeof( m ) ); }
    void init( ) { for( int i = 1; i <= sz; ++ i ) m[i][i] = 1; }

    friend Matrix operator * ( const Matrix &a, const Matrix &b ) {
        Matrix res( a.sz );
        for( int i = 1; i <= res.sz; ++ i ) {
            for( int k = 1; k <= res.sz; ++ k ) {
                if( a.m[i][k] == 0 ) continue;
                for( int j = 1; j <= res.sz; ++ j )
                    res.m[i][j] = ( res.m[i][j] + a.m[i][k] * b.m[k][j] ) % mod;
            }
        }
        return res;
    }
};

Matrix mqpow( Matrix a, ll k ) 
{
    Matrix res( a.sz ); res.init( );
    for( ; k; k >>= 1, a = a * a )
        if( k & 1 ) res = res * a;
    return res;
}

int main( )
{
    ios::sync_with_stdio( 0 ); cin.tie( 0 );
    
    cin >> p >> q >> a1 >> a2 >> n >> mod;

    if( n == 1 ) { cout << a1 << '\n'; return 0; }
    if( n == 2 ) { cout << a2 << '\n'; return 0; }
    
    Matrix T( 2 );
    T.m[1][1] = p, T.m[1][2] = q, T.m[2][1] = 1;
    T = mqpow( T, n - 2 );

    cout << ( T.m[1][1] * a2 + T.m[1][2] * a1 ) % mod << '\n';
    return 0;
}

Ex.2 P3216 数学作业

1. 题面梗概

将 $1 \dots n$ 顺序拼接数，求该数对 $m$ 取模的结果。

$1 \le n \le 10^{18}, 1 \le m \le 10^9$

2. 逻辑推导

我们可以简单的用 DP 思路写出转移方程。

设 $f_{i}$ 表示前 $i$ 个数拼接的结果，当前加入的数 $i$ 有 $d$ 位（ $d = \lfloor \lg i \rfloor + 1$ ），那么：

$$f_{i} = f_{i-1} \cdot 10^{d} + i \pmod m$$

但是问题在于 $n \leq 10^{18}$，这种线性的递推是 $O(n)$ 的，直接 TLE 一波送走。

这同时也是一个极具线性关系的递推式，所以我们可以想到使用矩阵快速幂加速递推。

首先构造出答案矩阵（这里是向量）：$\begin{bmatrix}dp[i] \\i \\1 \\ \end{bmatrix}$，接着用待定系数法构造出转移矩阵：$\begin{bmatrix}10^d & 1 & 0 \\0 & 1 & 1 \\0 & 0 & 1 \\\end{bmatrix}$，然后递推即可。

答案矩阵是怎么构造的呢？

答案矩阵一般都是单列的：

1. 先把 $dp[i]$ 写上，方便在最后直接输出。

$$\begin{bmatrix} dp[i] \\ \\ \\ \end{bmatrix}$$

2. 回顾我们的递推方程，$dp[i]$ 需要通过 $dp[i-1]$ 和 $i$ 推出，这两项都需要体现在矩阵里，为了转移，我们还需要 $i$。

$$\begin{bmatrix} dp[i] \\ i \\ \\ \end{bmatrix}$$

3. 由 $i-1$ 是如何推出 $i$ 的呢？所以矩阵里还需要一个常数 $1$。

小问题

我们发现 $d$ 会随着 $i$ 改变，这在矩阵快速幂中是不被允许的，因此不能一次矩阵快速幂解决，我们需要按照 $i$ 的位数，分别对每一段进行转移。

3. 代码实现

▶

AcCode

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

const int maxsz = 5;
const int maxn = 100 + 5;

ll dp[maxn];
ll n, mod;

struct Matrix {
    ll m[maxsz][maxsz];
    int sz;
    Matrix( int s = 0 ) : sz( s ) { memset( m, 0, sizeof( m ) ); }
    void init( ) { for( int i = 1; i <= sz; ++ i ) m[i][i] = 1; }

    friend Matrix operator * ( const Matrix &a, const Matrix &b ) {
        Matrix res( a.sz );
        for( int i = 1; i <= res.sz; ++ i ) {
            for( int k = 1; k <= res.sz; ++ k ) {
                if( a.m[i][k] == 0 ) continue;
                for( int j = 1; j <= res.sz; ++ j )
                    res.m[i][j] = ( res.m[i][j] + a.m[i][k] * b.m[k][j] ) % mod;
            }
        }
        return res;
    }
};

Matrix mqpow( Matrix a, ll k ) 
{
    Matrix res( a.sz ); res.init( );
    for( ; k; k >>= 1, a = a * a )
        if( k & 1 ) res = res * a;
    return res;
}

ll qpow( ll a, ll k ) {
    ll res = 1;
    for( ; k; k >>= 1, a = a * a % mod )
        if( k & 1 ) res = res * a % mod;
    return res;
}

int main( )
{
    ios::sync_with_stdio( 0 ); cin.tie( 0 );
    
    cin >> n >> mod;

    ll ans[4] = { 0, 0, 0, 1 };
    
    for( ll d = 1, st = 1; st <= n; d ++, st *= 10 )
    {
        ll ed;
        if( st > n / 10 ) ed = n;
        else ed = st * 10 - 1;

        ll cnt = ed - st + 1;

        Matrix T( 3 );
        T.m[1][1] = qpow( 10, d );
        T.m[1][2] = 1; T.m[1][3] = 1;
        T.m[2][1] = 0; T.m[2][2] = 1; T.m[2][3] = 1;
        T.m[3][1] = 0; T.m[3][2] = 0; T.m[3][3] = 1;
        T = mqpow( T, cnt );
        
        ll nxt[4] = { 0, 0, 0, 0 };
        for( int i = 1; i <= 3; ++ i )
            for( int j = 1; j <= 3; ++ j )
                nxt[i] = ( nxt[i] + ans[j] * T.m[i][j] ) % mod;
        
        for( int i = 1; i <= 3; ++ i ) ans[i] = nxt[i];
        if( st > n / 10 ) break;
    }

    cout << ans[1] % mod << '\n';
    return 0;
}

二、图上转移

Ex.1 P4159 迷路

1. 题面梗概

给定一张 $n$ 个点的有向图，边权为 $1 \sim 9$ 之间的整数。求从 1 号点出发，恰好在 $t$ 时刻到达 $n$ 号点的方案数。

$n \le 10, t \le 10^9, w_{i,j} \in [0, 9]$

2. 逻辑推导

在解决这道题之前，我们需要理解 “单位权图路径计数” 的本质。

首先我们先用计数 $DP$ 来演示：$dp[i][j] = \sum dp[i][k] + dp[k][j]$。

那么矩阵呢？设邻接矩阵为 $A$，其中 $A_{i.j} = 1$ 表示存在 $i \to j$ 的边，当我们计算 $A^2$ 时，根据乘法定义：

$$(A^2)_{i, j} = \sum_{k=1}^{n} A_{i, k} \cdot A_{k, j}$$

也就是说：

• $A_{i,k} \cdot A_{k,j} = 1$ 当且仅当存在路径 $i \to k$ 且存在路径 $k \to j$。

• 这意味着存在一条经过中间点 $k$、长度为 2 的路径 $i \to k \to j$。

• 对所有可能的中间点 $k$ 求和，得到的就是从 $i$ 到 $j$ 长度恰好为 2 的路径总数。

以此类推，$(A^t)_{i,j}$ 的意义就是：从 $i$ 到 $j$ 恰好经过 $t$ 条单位边（耗时为 $t$ ）的路径方案数。

在这里，矩阵乘法本质上是一种并行化的 DP。

特殊处理：边权不为 1

回到这一题，发现边权 $w \in [1, 9]$，意味着走过一条边需要耗费多个时间单位，没法直接用矩阵幂乘。

但是，由于边权 $w$ 极小，我们可以考虑类似 “分层图” 的技巧，取把原始节点拆成多个状态节点。

对于原图中的每一个点 $u$，我们将其拆分为 $9$ 个虚拟节点：

$$(u, 0), (u, 1), \dots, (u, 8)$$

定义：$(u, k)$ 表示 “处于通往 $u$ 的路径上，且还需要 $k$ 秒才能到达点 $u$ ”。（ $(u, 0)$ 表示已经到这个点了 ）

内部转移：在同一个点的虚拟节点间连边 $(u, k) \to (u, k-1)$，权值为 $1$。

外部转移：对于边 $(u, v, w)$，我们连边 $(u, 0) \to (v, w - 1)$，权值为

3. 代码实现

▶

AcCode

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

const int maxsz = 100 + 5;
const int maxn = 100 + 5;
const int mod = 2009;

int n, t;

struct Matrix {
    ll m[maxsz][maxsz];
    int sz;
    Matrix( int s = 0 ) : sz( s ) { memset( m, 0, sizeof( m ) ); }
    void init( ) { for( int i = 1; i <= sz; ++ i ) m[i][i] = 1; }

    friend Matrix operator * ( const Matrix &a, const Matrix &b ) {
        Matrix res( a.sz );
        for( int i = 1; i <= res.sz; ++ i ) {
            for( int k = 1; k <= res.sz; ++ k ) {
                if( a.m[i][k] == 0 ) continue;
                for( int j = 1; j <= res.sz; ++ j )
                    res.m[i][j] = ( res.m[i][j] + a.m[i][k] * b.m[k][j] ) % mod;
            }
        }
        return res;
    }
};

Matrix mqpow( Matrix a, ll k ) 
{
    Matrix res( a.sz ); res.init( );
    for( ; k; k >>= 1, a = a * a )
        if( k & 1 ) res = res * a;
    return res;
}

char c;

int getid( int i, int j ) { return i + j * n; }

int main( )
{
    ios::sync_with_stdio( 0 ); cin.tie( 0 );

    cin >> n >> t;
    int sz = n * 9;
    Matrix A( sz );

    for( int i = 1; i <= n; ++ i )
    {
        for( int j = 1; j <= 8; ++ j )
            A.m[getid( i, j )][getid( i, j - 1 )] = 1;
        
        for( int j = 1; j <= n; ++ j )
        {
            cin >> c;
            if( c == '0' ) continue;
            int w = c - '0';
            A.m[getid( i, 0 )][getid( j, w - 1 )] = 1;
        }
    }

    Matrix res = mqpow( A, t );
    cout << res.m[getid( 1, 0 )][getid( n, 0 )] << '\n';
    return 0;
}

Ex.2 P2151 HH 去散步

1. 题面梗概

给定一个有向图，求从 $A$ 到 $B$ 长度为 $l$ 的路径数量。

限制条件：不能沿着刚刚来过来的路走回去（不能走回头路，仅限 1 步）。

$N \leq 50, M \leq 60, t \leq 2^{30}$，结果模 $45989$。

2. 逻辑推导

如果只是求从 $A$ 到 $B$ 的路径数量，通过对邻接矩阵的幂乘可以直接通过 $O( \log N)$ 的复杂度快速求得。

但是这题有一个限制，不能立即原路返回，意味着当前的决策依赖于 “上一步是从那条边过来的”，在传统的以点为状态的模型中，这就产生了 后效性，没法通过矩阵乘法实现。

PS：矩阵乘法本质是高并发的 DP，因此必须具备无后效性。

点边转换

注意到我们要记录上一条边，不妨把状态记录从 点，转化为 边，记录刚刚走过了哪条有向边。

那么我们可以把原图上的 $M$ 条无向边，拆解为 $2M$ 条有向边。

这里有一个小技巧，Cue 链式前向星，i ^ 1 可以快速得到对应的反向边，前提是 0 - index，具体推导就略了。

接着定义状态：$dp[t][i]$ 表示第 $t$ 步走过有向边 $i$ 的方案数。

转移矩阵构造

对于两条有向边 $X(u \to v)$，$Y(p \to q)$，能够从 $X$ 转移到 $Y$ 当且仅当：

1. 边 $X$ 的终点等于边 $Y$ 的起点。
2. 边 $Y$ 不是边 $X$ 的反向边。

满足上述条件时，转移矩阵 $T[i][j] = 1$。这是个 $2M \times 2M$ 的矩阵，最大规模 $120$，复杂度还是很简单的。

初始状态

由于第一步没有前驱状态，不受限制，所以我们设置初始向量 $V$ 记录所有以起点 $A$ 为起点的有向边的方案数为 $1$。

然后，总步数因为第一步已经走完，所以矩阵加速的步数为 $t - 1$。

最后统计答案：终点为 $B$ 的有向边的方案数。

3. 代码实现

▶

AcCode

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

const int maxm = 60 + 5;
const int maxsz = 2 * maxm;
const int mod = 45989;

struct Edge { int u, v; } e[maxm << 1];
int cnt = 0;

int n, m, t, st, ed;

struct Matrix {
    ll m[maxsz][maxsz];
    int sz;
    Matrix( int s = 0 ) : sz( s ) { memset( m, 0, sizeof( m ) ); }
    void init( ) { for( int i = 0; i < sz; ++ i ) m[i][i] = 1; }

    friend Matrix operator * ( const Matrix &a, const Matrix &b ) {
        Matrix res( a.sz );
        for( int i = 0; i < res.sz; ++ i ) {
            for( int k = 0; k < res.sz; ++ k ) {
                if( a.m[i][k] == 0 ) continue;
                for( int j = 0; j < res.sz; ++ j )
                    res.m[i][j] = ( res.m[i][j] + a.m[i][k] * b.m[k][j] ) % mod;
            }
        }
        return res;
    }
};

Matrix mqpow( Matrix a, ll k ) 
{
    Matrix res( a.sz ); res.init( );
    for( ; k; k >>= 1, a = a * a )
        if( k & 1 ) res = res * a;
    return res;
}

int main( )
{
    ios::sync_with_stdio( 0 ); cin.tie( 0 );
    
    cin >> n >> m >> t >> st >> ed;
    for( int i = 0; i < m; ++ i )
    {
        int u, v; cin >> u >> v;
        e[cnt ++] = { u, v };
        e[cnt ++] = { v, u };
    }

    Matrix T( cnt );

    for( int i = 0; i < cnt; ++ i )
    {
        for( int j = 0; j < cnt; ++ j )
        {
            if( e[i].v == e[j].u && j != (i ^ 1) )
                T.m[i][j] = 1;
        }
    }

    Matrix V( cnt );
    for( int i = 0; i < cnt; ++ i )
    {
        if( e[i].u == st ) V.m[0][i] = 1;
    }

    Matrix res = V * mqpow( T, t - 1 );

    ll ans = 0;
    for( int i = 0; i < cnt; ++ i )
    {
        if( e[i].v == ed ) ans = ( ans + res.m[0][i] ) % mod;
    }

    cout << ans << '\n';
    return 0;
}

Ex.3 ABC445F - Exactly K Steps 2

1. 题面梗概

给一个 $N$ 个点的完全有向图，每条边 $(i, j)$ 都有代价 $C_{i,j}$。给定 $K$，对于每个起点 $s$，求出恰好经过 $K$ 条边并回到 $s$ 的最小总代价。

$N \leq 100, K \leq 1e9$

2. 逻辑推导

考虑 DP，设 $dp[k][i][j]$ 表示恰好经过 $k$ 条边从 $i$ 到 $j$ 的最小代价。

那么很容易能写出一个转移方程。

$$dp[k+1][i][j] = \min_{1 \le m \le N} \{ dp[k][i][m] + C_{m,j} \}$$

为了走 $k + 1$ 步到 $j$，我可以先走 $k$ 步到某个中间点 $m$，然后再花 $C_{m,j}$ 的代价走最后一步。

但是！$k \leq 10^9$，如果直接线性状态转移，绝对会 TLE ，我们需要优化。

对于这个数量级的规模，最简单想到的一种优化就是 二分 或者 倍增 这种带 $\log$ 的复杂度。

我们再拆解一下这个转移方程：

$$dp[a + b][i][j] = \min_{1 \leq p \leq n}(dp[a][i][p] + dp[b][p][j])$$

有思路了吗？没有？我们再对比一下下面这个柿子：

$$C[i][j] = \sum( A[i][p] \cdot B[p][j])$$

是不是惊人的相似，这里解释一下，下面这个是 矩阵乘法，而前面那个状态转移其实是满足 结合律 的。

我们只需要把 求和 替换为 取最小值，乘法 换成 加法，这个问题就可以用矩阵快速幂优化。

插入一句，这玩意学名好像叫 $(min, +)$ 卷积。

快速幂

既然这种运算符合 结合律 的，我们就可以定义一个初始代价矩阵 $M$，其中 $M_{i,j} = C_{i, j}$，经过 $K$ 条边的最短路径矩阵就是 $M^{K}$。

单位矩阵：在这个 $(min, +)$ 运算下，单位矩阵 $I$ 满足 $I_{i,i} = 0$，其余元素为 $\infty$。

广义单位矩阵

这里插入一节拓展知识。

在标准的算术矩阵中，单位矩阵 $I$ 的对角线是 $1$（乘法单位元），其余是 $0$（加法单位元）。

但在 $(\min, +)$ 代数系中：

• 乘法 $\to$ 加法：其单位元是 $0$（因为 $x + 0 = x$）。

• 加法 $\to$ 最小值：其单位元是 $\infty$（因为 $\min(x, \infty) = x$）。

因此，$(\min, +)$ 意义下的单位矩阵 $I$ 必须满足：

• 对角线 $I_{i,i} = 0$：代表走 $0$ 步的代价为 $0$。

• 其余 $I_{i,j} = \infty$：代表走 $0$ 步无法到达其他位置。

但是这里我偷懒了，我直接把矩阵初始化为第一步的情况，就避免了 走 0 步 的边界处理。

3. 代码实现

▶

AcCode

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

const int maxn = 100 + 5;
const ll inf = 1e18;

ll n, k;

struct Mat {
    ll mat[maxn][maxn];
    Mat( ) { memset( mat, 0x3f, sizeof (mat) ); }
    
    // (min, +) 矩阵乘法
    Mat mul( const Mat &A, const Mat &B )
    {
        Mat res;
        for( int k = 1; k <= n; ++ k )
        {
            for( int i = 1; i <= n; ++ i )
            {
                if( A.mat[i][k] == inf ) continue;
                for( int j = 1; j <= n; ++ j )
                {
                    if( B.mat[k][j] == inf ) continue;
                    res.mat[i][j] = min( res.mat[i][j], A.mat[i][k] + B.mat[k][j] );
                }
            }
        }
        return res;
    }
    // 矩阵快速幂
    Mat pow( Mat A, ll k )
    {
        Mat res = A;
        while( k )
        {
            if( k & 1 ) res = mul( res, A );
            A = mul( A, A );
            k >>= 1;
        }
        return res;
    }
};

void solve( )
{
    cin >> n >> k;

    Mat c;
    for( int i = 1; i <= n; ++ i )
    {
        for( int j = 1; j <= n; ++ j )
        {
            cin >> c.mat[i][j];
        }
    }

	// 刚开始是 1 步，还需要再走 k - 1 步。
    Mat ans = c.pow( c, k - 1 );

    for( int i = 1; i <= n; ++ i )
    {
        cout << ans.mat[i][i] << '\n';
    }
}

int main( )
{
    ios::sync_with_stdio(0);
    cin.tie(0); cout.tie(0);
    solve( );
    return 0;
}