矩阵快速幂优化 DP（最全！）

156****1682

2026-07-19 21:30:27

发布于：北京

11阅读

0回复

0点赞

矩阵快速幂优化 DP

1. 矩阵快速幂优化 DP 的定义

若目标 DP 满足线性和马尔可夫性，存在 $k$ 维状态向量 $V_i$ 和 $k \times k$ 转移矩阵 $M$ ，满足

$V_{i+1} = M \cdot V_i$

则

$V_n = M^{\,n-k} \cdot V_k$

因为矩阵乘法具有结合律，所以可以利用快速幂将 DP 的时间复杂度优化至 $O(k^3 \log n)$ 。

2. 一般线性递推的矩阵构建

对于递推

$dp[i] = c_1 dp[i-1] + c_2 dp[i-2] + \dots + c_k dp[i-k]$

状态向量为

$V_i = [dp[i],\ dp[i-1],\ \dots,\ dp[i-k+1]]^{\!T}$

则转移矩阵为

$M = \begin{bmatrix} c_1 & c_2 & \dots & c_{k-1} & c_k \\ 1 & 0 & \dots & 0 & 0 \\ 0 & 1 & \dots & 0 & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & \dots & 1 & 0 \end{bmatrix}_{k \times k}$

初始向量取 $V_{k-1} = [dp[k-1],\ dp[k-2],\ \dots,\ dp[0]]^{T}$ 。则

$dp[n] = \bigl( M^{\,n-k+1} \cdot V_{k-1} \bigr)[0]$

示例代码:

#include<iostream>
using namespace std ;
const int K = 2 ;
const long long MOD = 1e9+7 ;
// c = a * b，并取模
void mul(long long a[K][K], long long b[K][K], long long c[K][K]){
    for(int i = 0 ; i < K ; i++){
        for(int j = 0 ; j < K ; j++){
            long long sum = 0 ;
            for(int t = 0 ; t < K ; t++){
                sum = (sum + a[i][t] * b[t][j]) % MOD ;
            }
            c[i][j] = sum ;
        }
    }
}
// ans = M^n
void pow(long long M[K][K] , long long n , long long ans[K][K]) {
    long long p[K][K] , t[K][K] ;
    for (int i = 0 ; i < K ; i++){
        for (int j = 0 ; j < K ; j++){
            p[i][j] = M[i][j] ;
        }
    }
    // 初始化 ans 为单位矩阵
    for(int i = 0 ; i < K ; i++){
        for(int j = 0 ; j < K ; j++){
            ans[i][j] = (i == j) ;
        }
    }
    while(n){
        if(n & 1){
            mul(ans, p, t) ;
            for(int i = 0 ; i < K ; i++){
                for(int j = 0 ; j < K ; j++){
                    ans[i][j] = t[i][j] ;
                }
            }
        }
        n >>= 1 ;
        mul(p, p , t) ;
        for(int i = 0 ; i < K ; i++){
            for(int j = 0 ; j < K ; j++){
                p[i][j] = t[i][j] ;
            }
        }
    }
}

int main() {
    long long n ;
    cin >> n ;
    if(n <= 1){
        cout << n << endl ;
        return 0 ;
    }
    // 以斐波那契数列为例
    // dp[i] = 1*dp[i-1] + 1*dp[i-2]
    // 转移矩阵: M = 1 1
    //              1 0
    long long M[K][K] = {{1 , 1} , {1 , 0}} ;
    long long res[K][K] ;
    pow(M , n - 1 , res) ;
    // 斐波那契数列的初始向量使得结果恰好为 res[0][0]
    // cout << res[0][0] << endl ;
    // 作为一般性示例，依然保留初始向量部分
    long long v[K] = {1 , 0} ;
    long long ans[K] = {0 , 0} ;
    // 乘初始向量
    for(int i = 0 ; i < K ; i++){
        for(int j = 0 ; j < K ; j++){
            ans[i] = (ans[i] + res[i][j] * v[j]) % MOD ;
        }
    }
    cout << ans[0] << endl;
    return 0;
}

3. 非齐次多项式驱动项的增广矩阵

设递推包含 $d$ 次多项式驱动项 $P(i)$ ：

$dp[i] = \sum_{j=1}^k c_j\, dp[i-j] + P(i)$

将 $P(i+1)$ 按 $i$ 的升幂展开为 $P(i+1)=q_0+q_1 i+q_2 i^2+\dots+q_d i^d$ 。取多项式基底为升幂排列，构造增广状态向量

$V_i = [dp[i],\ dp[i-1],\ \dots,\ dp[i-k+1],\ 1,\ i,\ i^2,\ \dots,\ i^d]^{T}$

总维度为 $k+d+1$ 。递推写作 $V_{i+1}=M\cdot V_i$ ，转移矩阵具有分块形式

$M = \left[ \begin{array}{cccc|cccc} c_1 & c_2 & \cdots & c_k & q_0 & q_1 & q_2 & \cdots & q_d \\ 1 & 0 & \cdots & 0 & 0 & 0 & 0 & \cdots & 0 \\ 0 & 1 & \cdots & 0 & 0 & 0 & 0 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots & \vdots & \vdots & & \vdots \\ 0 & 0 & \cdots & 1 & 0 & 0 & 0 & \cdots & 0 \\ \hline 0 & 0 & \cdots & 0 & \binom{0}{0} & 0 & 0 & \cdots & 0 \\ 0 & 0 & \cdots & 0 & \binom{1}{0} & \binom{1}{1} & 0 & \cdots & 0 \\ 0 & 0 & \cdots & 0 & \binom{2}{0} & \binom{2}{1} & \binom{2}{2} & \cdots & 0 \\ \vdots & \vdots & & \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & 0 & \binom{d}{0} & \binom{d}{1} & \binom{d}{2} & \cdots & \binom{d}{d} \end{array} \right]_{(k+d+1)\times(k+d+1)}$

右下角块为 $(d+1)\times(d+1)$ 下三角矩阵，其中的第 $p$ 行第 $t$ 列为 $\binom{p}{t}$ ， $p,t=0,1,\dots,d$ ，当 $t>p$ 时元素为 $0$ 。该块将 $[1,i,i^2,\dots,i^d]^T$ 映射为 $[1,i+1,(i+1)^2,\dots,(i+1)^d]^T$ 。

递推从 $i=1$ 开始，所有负数下标的 $dp[\cdot]$ 均视为 $0$ 。初始向量取

$V_1 = [dp[1],\ dp[0],\ 0,\ \dots,\ 0,\ 1,\ 1,\ 1^2,\ \dots,\ 1^d]^{T}$

则

$dp[n] = \bigl( M^{\,n-1} \cdot V_1 \bigr)[0]$

示例代码：

#include<iostream>
using namespace std ;
const int K = 5 ; // k = 2 , d = 2  ->  K = k + d + 1 = 5
const long long MOD = 1e9+7 ;
// c = a * b，并取模
void mul(long long a[K][K], long long b[K][K], long long c[K][K]){
    for(int i = 0 ; i < K ; i++){
        for(int j = 0 ; j < K ; j++){
            long long sum = 0 ;
            for(int t = 0 ; t < K ; t++){
                sum = (sum + a[i][t] * b[t][j]) % MOD ;
            }
            c[i][j] = sum ;
        }
    }
}
// ans = M^n
void pow(long long M[K][K] , long long n , long long ans[K][K]) {
    long long p[K][K] , t[K][K] ;
    for (int i = 0 ; i < K ; i++){
        for (int j = 0 ; j < K ; j++){
            p[i][j] = M[i][j] ;
        }
    }
    // 初始化 ans 为单位矩阵
    for(int i = 0 ; i < K ; i++){
        for(int j = 0 ; j < K ; j++){
            ans[i][j] = (i == j) ;
        }
    }
    while(n){
        if(n & 1){
            mul(ans, p, t) ;
            for(int i = 0 ; i < K ; i++){
                for(int j = 0 ; j < K ; j++){
                    ans[i][j] = t[i][j] ;
                }
            }
        }
        n >>= 1 ;
        mul(p, p , t) ;
        for(int i = 0 ; i < K ; i++){
            for(int j = 0 ; j < K ; j++){
                p[i][j] = t[i][j] ;
            }
        }
    }
}

int main() {
    // 递推式: dp[i] = 2dp[i-1] + dp[i-2] + 3i^2
    // 设 dp[0] = 0 , dp[1] = 1 
    long long n ;
    cin >> n ;
    if(n == 0){cout << 0 << endl ;return 0 ;}
    if(n == 1){cout << 1 << endl ;return 0 ;}
    // 状态向量 V_i = [dp[i], dp[i-1], 1, i, i^2]^T
    // P(i+1) = 3(i+1)^2 = 3 + 6i + 3i^2 -> q0=3, q1=6, q2=3
    // 转移矩阵，按照第 3 节构造
    long long M[K][K] = {
        {2, 1, 3, 6, 3},   // 新 dp
        {1, 0, 0, 0, 0},   // dp 下移
        {0, 0, 1, 0, 0},   // 1 不变
        {0, 0, 1, 1, 0},   // i -> i+1 = 1+i
        {0, 0, 1, 2, 1}    // i^2 -> (i+1)^2 = 1+2i+i^2
    };
    long long res[K][K] ;
    pow(M, n - 1, res) ;
    // 初始向量 V_1 (i=1): dp[1]=1, dp[0]=0, 1, i=1, i^2=1
    long long v[K] = {1, 0, 1, 1, 1} ;
    long long ans[K] = {0} ;
    // 乘初始向量
    for(int i = 0 ; i < K; i++){
        for(int j = 0 ; j < K; j++){
            ans[i] = (ans[i] + res[i][j] * v[j]) % MOD ;
        }
    }
    cout << ans[0] << endl ;
    return 0;
}

4. 前缀和累加器的嵌入矩阵

对线性递推 $dp[i] = \sum_{j=1}^k c_j\, dp[i-j]$ ，需要同时维护前缀和 $S[i] = \sum_{t=0}^i dp[t]$ 。定义增广状态向量

$V_i = [dp[i],\ dp[i-1],\ \dots,\ dp[i-k+1],\ S[i]]^{T}$

维度为 $k+1$ 。递推关系对 $i \ge k$ 成立，转移矩阵为

$M = \begin{bmatrix} c_1 & c_2 & \dots & c_k & 0 \\ 1 & 0 & \dots & 0 & 0 \\ 0 & 1 & \dots & 0 & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & \dots & 1 & 0 \\ \hline c_1 & c_2 & \dots & c_k & 1 \end{bmatrix}_{(k+1)\times(k+1)}$

最后一行复制系数 $c_1,\dots,c_k$ 以计算新的 $dp[i+1]$ ，末列的 $1$ 将旧前缀和 $S[i]$ 累加到 $S[i+1]$ 中。

初始向量取

$V_{k-1} = [dp[k-1],\ dp[k-2],\ \dots,\ dp[0],\ S[k-1]]^{T}$

其中 $S[k-1] = \sum_{t=0}^{k-1} dp[t]$ 。则对 $n \ge k-1$ ，

$S[n] = \bigl( M^{\,n-k+1} \cdot V_{k-1} \bigr)[k]$

示例代码：

#include<iostream>
using namespace std ;
const int K = 3 ; // k = 2 + 前缀和嵌入 -> K = k + 1 = 3
const long long MOD = 1e9+7 ;
// c = a * b，并取模
void mul(long long a[K][K], long long b[K][K], long long c[K][K]){
    for(int i = 0 ; i < K ; i++){
        for(int j = 0 ; j < K ; j++){
            long long sum = 0 ;
            for(int t = 0 ; t < K ; t++){
                sum = (sum + a[i][t] * b[t][j]) % MOD ;
            }
            c[i][j] = sum ;
        }
    }
}

// ans = M^n
void pow(long long M[K][K] , long long n , long long ans[K][K]) {
    long long p[K][K] , t[K][K] ;
    for (int i = 0 ; i < K ; i++){
        for (int j = 0 ; j < K ; j++){
            p[i][j] = M[i][j] ;
        }
    }
    // 初始化 ans 为单位矩阵
    for(int i = 0 ; i < K ; i++){
        for(int j = 0 ; j < K ; j++){
            ans[i][j] = (i == j) ;
        }
    }
    while(n){
        if(n & 1){
            mul(ans, p, t) ;
            for(int i = 0 ; i < K ; i++){
                for(int j = 0 ; j < K ; j++){
                    ans[i][j] = t[i][j] ;
                }
            }
        }
        n >>= 1 ;
        mul(p, p , t) ;
        for(int i = 0 ; i < K ; i++){
            for(int j = 0 ; j < K ; j++){
                p[i][j] = t[i][j] ;
            }
        }
    }
}

int main() {
    // 斐波那契数列求和
    // dp[i] = dp[i-1] + dp[i-2]
    // dp[0] = 0 , dp[1] = 1
    // 状态向量 V_i = [dp[i], dp[i-1], S[i]]^T
    long long n ;
    cin >> n ;
    if(n == 0){cout << 0 << endl ;return 0 ;}
    if(n == 1){cout << 1 << endl ;return 0 ;}
    // 转移矩阵，按照第 4 节构造，k=2, c1=1, c2=1
    long long M[K][K] = {
        {1, 1, 0},   // 新 dp[i+1] = 1*dp[i] + 1*dp[i-1]
        {1, 0, 0},   // dp[i] 下移
        {1, 1, 1}    // 新 S[i+1] = 新dp + S[i] = 1*dp[i]+1*dp[i-1] + 1*S[i]
    };
    long long res[K][K] ;
    pow(M, n - 1, res) ;
    // 初始向量 V_1 (i=1): dp[1]=1, dp[0]=0, S[1]=dp[1]=1
    long long v[K] = {1, 0, 1} ;
    long long ans[K] = {0} ;
    // 乘初始向量
    for(int i = 0 ; i < K; i++){
        for(int j = 0 ; j < K; j++){
            ans[i] = (ans[i] + res[i][j] * v[j]) % MOD ;
        }
    }
    // 前缀和存于 ans[2]
    cout << ans[2] << endl ;
    return 0;
}

5. 高阶依赖与拆点构造

当问题依赖过去的状态，仅记录当前点将无法判断下一步的合法性，需要将状态重新定义为边，以补齐所依赖的历史信息，恢复马尔可夫性。在此类图问题中，前几步所走的边是所需历史信息的载体。

设边集为 $E$ ，状态向量

$V_i = [\mathrm{cnt}(e_1),\ \mathrm{cnt}(e_2),\ \dots,\ \mathrm{cnt}(e_E)]^{T}$

其中 $\mathrm{cnt}(e)$ 表示第 $i$ 步结束时恰落在边 $e$ 终点的路径数。转移矩阵 $M$ 是 $E\times E$ 的 $0\text{-}1$ 矩阵，其元素定义为

$M_{e',\,e} = \begin{cases} 1, & \text{若 } \mathrm{end}(e) = \mathrm{start}(e') \text{ 且满足题设附加约束} \\ 0, & \text{否则} \end{cases}$

这里 $\mathrm{start}(e),\mathrm{end}(e)$ 分别表示边 $e$ 的起点与终点。递推式为 $V_{i+1} = M \cdot V_i$ 。

初始向量 $V_1$ 中，以起点出发的边对应的分量为 $1$ ，其余分量为 $0$ 。则长度为 $n$ 的路径中，到达目标顶点 $t$ 的路径总数为

$\mathrm{ans} = \sum_{e:\,\mathrm{end}(e)=t} \bigl( M^{\,n-1} \cdot V_1 \bigr)_e$

#include<iostream>
using namespace std ;
const int K = 4 ; // 边的数量 E = 4
const long long MOD = 1e9+7 ;
// c = a * b，并取模
void mul(long long a[K][K], long long b[K][K], long long c[K][K]){
    for(int i = 0 ; i < K ; i++){
        for(int j = 0 ; j < K ; j++){
            long long sum = 0 ;
            for(int t = 0 ; t < K ; t++){
                sum = (sum + a[i][t] * b[t][j]) % MOD ;
            }
            c[i][j] = sum ;
        }
    }
}
// ans = M^n
void pow(long long M[K][K] , long long n , long long ans[K][K]) {
    long long p[K][K] , t[K][K] ;
    for (int i = 0 ; i < K ; i++){
        for (int j = 0 ; j < K ; j++){
            p[i][j] = M[i][j] ;
        }
    }
    // 初始化 ans 为单位矩阵
    for(int i = 0 ; i < K ; i++){
        for(int j = 0 ; j < K ; j++){
            ans[i][j] = (i == j) ;
        }
    }
    while(n){
        if(n & 1){
            mul(ans, p, t) ;
            for(int i = 0 ; i < K ; i++){
                for(int j = 0 ; j < K ; j++){
                    ans[i][j] = t[i][j] ;
                }
            }
        }
        n >>= 1 ;
        mul(p, p , t) ;
        for(int i = 0 ; i < K ; i++){
            for(int j = 0 ; j < K ; j++){
                p[i][j] = t[i][j] ;
            }
        }
    }
}
int main() {
    // 示例问题：
    // e0: A->B (红)
    // e1: A->C (蓝)
    // e2: B->C (蓝)
    // e3: C->A (红)
    // 不能连续走同色的边
    // 求从 A 出发，走 n 步，回到 A 的路径数
    long long n ;
    cin >> n ;
    if(n == 0 || n == 1){cout << 0 << endl ;return 0 ;}
    // 转移矩阵 M[i][j] 表示是否可以从边 j 走到边 i，可以为1，不可以为0
    // 行i，列j
    long long M[K][K] = {
        {0, 0, 0, 0},  // e0
        {0, 0, 0, 1},  // e1
        {1, 0, 0, 0},  // e2
        {0, 1, 1, 0}   // e3
    };
    long long res[K][K] ;
    pow(M, n - 1, res) ;
    // 初始向量 V_1：第一步从 A 出发，可走 e0 或 e1
    long long v[K] = {1, 1, 0, 0} ;
    long long ans[K] = {0} ;
    // 乘初始向量
    for(int i = 0 ; i < K; i++){
        for(int j = 0 ; j < K; j++){
            ans[i] = (ans[i] + res[i][j] * v[j]) % MOD ;
        }
    }
    // 终点为 A 的边是 e3
    cout << ans[3] << endl ;
    return 0;
}

6. 多维状态压缩的向量化矩阵

若 DP 的第二维是有限状态 $mask \in \{0,1,\dots,L-1\}$ ，且转移形如

$dp[i][\mathrm{new}] = \sum_{\mathrm{old}} w(\mathrm{old} \to \mathrm{new})\; dp[i-1][\mathrm{old}]$

则可视为一阶递推，将第二维展平为向量

$V_i = [dp[i][0],\ dp[i][1],\ \dots,\ dp[i][L-1]]^{T}$

转移矩阵 $M$ 是 $L\times L$ 的邻接矩阵，其元素为

$M[\mathrm{new}][\mathrm{old}] = w(\mathrm{old} \to \mathrm{new})$

权值 $w$ 直接编码状态机的一切转移，常取 $0$ 或 $1$ 。递推写作 $V_i = M \cdot V_{i-1}$ 。

初始向量 $V_0$ 由初始分布 $dp[0][mask]$ 直接给出。则

$dp[n][mask] = \bigl( M^{\,n} \cdot V_0 \bigr)[mask]$

此式对应于一般形式中取 $k=1$ 的特例。

#include<iostream>
using namespace std ;
// 题目：
// 给定 N 个节点 M 条有向边（可能有重边），求从 S 到 T 恰好经过 k 条边的路径数，模 1e9+7。
// N <= 50, k <= 1e18
const int K = 55 ; // 最大节点数
const long long MOD = 1e9+7 ;
// c = a * b，并取模
void mul(long long a[K][K], long long b[K][K], long long c[K][K]){
    for(int i = 0 ; i < K ; i++){
        for(int j = 0 ; j < K ; j++){
            long long sum = 0 ;
            for(int t = 0 ; t < K ; t++){
                sum = (sum + a[i][t] * b[t][j]) % MOD ;
            }
            c[i][j] = sum ;
        }
    }
}
// ans = M^n
void pow(long long M[K][K] , long long n , long long ans[K][K]) {
    long long p[K][K] , t[K][K] ;
    for (int i = 0 ; i < K ; i++){
        for (int j = 0 ; j < K ; j++){
            p[i][j] = M[i][j] ;
        }
    }
    // 初始化 ans 为单位矩阵
    for(int i = 0 ; i < K ; i++){
        for(int j = 0 ; j < K ; j++){
            ans[i][j] = (i == j) ;
        }
    }
    while(n){
        if(n & 1){
            mul(ans, p, t) ;
            for(int i = 0 ; i < K ; i++){
                for(int j = 0 ; j < K ; j++){
                    ans[i][j] = t[i][j] ;
                }
            }
        }
        n >>= 1 ;
        mul(p, p , t) ;
        for(int i = 0 ; i < K ; i++){
            for(int j = 0 ; j < K ; j++){
                p[i][j] = t[i][j] ;
            }
        }
    }
}
int main(){
    int n , m ;
    long long k ;
    int s , t ;
    cin >> n >> m >> k >> s >> t ;
    // 构建转移矩阵 M (列 -> 行)
    // M[v][u] = 从节点 u 到节点 v 的路径数
    long long M[K][K] = {0} ;
    for(int i = 0 ; i < m ; i++){
        int u , v ;
        cin >> u >> v ;
        M[v][u] = M[v][u] + 1 ;
    }
    long long Mk[K][K] ;
    pow(M, k, Mk) ;
    // 初始向量 V0：只有起点 S 处为 1，其余为 0
    long long V0[K] = {0} ; V0[s] = 1 ;
    long long Vk[K] = {0} ;
    for(int i = 1 ; i <= n ; i++){
        for(int j = 1 ; j <= n ; j++){
            Vk[i] = (Vk[i] + Mk[i][j] * V0[j]) % MOD ;
        }
    }
    // 到达终点 T 的路径数
    cout << Vk[t] << endl ;
    return 0 ;
}