title: '[动态规划]矩阵加速'
date: 2026-03-31 21:33:43
tags: 动态规划 数论

矩阵加速DP转移

前置知识：矩阵快速幂

给定一个 $n$ 阶的矩阵 $A$ 以及一个非负整数 $k$ ，要计算矩阵 $A^{k}$ ，当 $k=0$ 时， $A^0$ 是 $n$ 阶单位矩阵 $I_n$ 。

即有 $A^{k} =A \times A \times ...\times A$ ，其中有 $k$ 个 $A$ 。

若 $k$ 数值小的时候，可以直接暴力计算。但是如果 $k\ge 10^{7}$ 的情况时，暴力做法必然会超时。

在此之前我们学过快速幂：即利用二进制的思想，将原来的底数不断扩大，做到 $O(log_2 k)$ 的时间复杂度求出一个数的 $k$ 次方。

i64 binpow(i64 a,i64 b,i64 m) {
  a%=m;
  i64 res=1;
  while(b>0) {
    if(b&1)res=res*a%m;
    a=a*a%m;
    b>>=1;
  }
  return res;
}

利用这个思想，我们可以把矩阵看作一个整体，套上整数快速幂的模板：

//创建Mat结构体，包含矩阵乘法
Mat binpow(Mat A,i64 x,i64 mod){
    Mat I;
    I.init_b();//初始化单位矩阵
    while(x){
        if(x&1) I=I*A;
        A=A*A;
        x>>=1;
    }
    return I;
}

求矩阵的 $k$ 次方这个问题就迎刃而解了。

矩阵加速

我们知道，动态规划可以看作一个递推的过程。其在遍历的过程中必然会推导出一个状态转移方程，只要整理清楚了这个状态转移方程，就可以对这个状态转移方程进行矩阵加速。举一个非常常见的例子：

斐波那契数列：

f(n) = \begin{cases} 0, & n = 0 \\ 1, & n = 1 \\ f(n-1) + f(n-2), & n \ge 2 \end{cases}

前面两个初始值不用太关心，我们主要关心当 $n \ge 2$ 时的情况：

f(n)=f(n-1)+f(n-2)

我们可以得到：

f(n+1)=f(n)+f(n-1)

假设向量 :

{I_n}^T=\begin{bmatrix} f(n) \\ f(n-1) \end{bmatrix}

{I_{n-1}}^T=\begin{bmatrix} f(n-1)\\ f(n-2) \end{bmatrix}

因为具有递推关系故存在一个 $2$ 阶矩阵 $A$ 于 ${I_{n-1}}^T$ 前面才能进行状态转移，即：

{I_n}^T=A \times {I_{n-1}}^T

也就是：

\begin{bmatrix} f(n) \\ f(n-1) \end{bmatrix}= \begin{bmatrix} a & b\\ c & d \end{bmatrix} \times \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}}^{k-1} \times \begin{bmatrix} f(2)\\ f(1) \end{bmatrix}

假设矩阵 $A$ 经过 $k-1$ 自乘之后变成了:

\begin{bmatrix} x & y\\ z & w \end{bmatrix}

由于矩阵乘法就有 $f(n)=x \times f(2)+y \times f(1)$ ;

矩阵自乘的过程可以运用矩阵快速幂。因此上式所求出的值即是第 $n$ 项斐波那契数列的值(一般会进行取模操作)。

这样一来我们就可以以 $O(log_2 n \times k^2)$ 的时间复杂度求出第 $n$ 项的各个值了。 $k$ 是矩阵的阶数。

例题

https://www.luogu.com.cn/problem/P3216

P3216 [HNOI2011] 数学作业

题目描述

小 C 数学成绩优异，于是老师给小 C 留了一道非常难的数学作业题：

给定正整数 $n,m$ ，要求计算 $\text{Concatenate}(n) \bmod \ m$ 的值，其中 $\text{Concatenate}(n)$ 是将 $1 \sim n$ 所有正整数顺序连接起来得到的数。

例如， $n = 13$ ， $\text{Concatenate}(n) = 12345678910111213$ 。小 C 想了大半天终于意识到这是一道不可能手算出来的题目，于是他只好向你求助，希望你能编写一个程序帮他解决这个问题。

【数据范围】

对于 $30\%$ 的数据， $1\le n \le 10^6$ ；
对于 $100\%$ 的数据， $1\le n \le 10^{18}$ ， $1\le m \le 10^9$ 。

不妨枚举一下：

$dp[1]=1,dp[2]=12$ ,即有 $dp[2]=dp[1] \times 10^{1}+2$ ，

$dp[2]=12,dp[3]=123$ ,即有 $dp[3]=dp[2]\times10^{1}+3，···$

$dp[9]=123456789,dp[10]=12345678910$ ,

有 $dp[10]=dp[9]*10^{2}+10$ ,

我们可以鲁莽的得到这个题目的 $DP$ 状态转移方程：

dp[i]=10^{len(i)}\times dp[i-1]+i

得出 $DP$ 状态转移方程就很简单了。但是题目中的 $1\le n \le 10^{18}$ ,数据特别大，可以考虑矩阵乘法。

对于 $len(i)$ ，在某一段范围是固定的，比如说， $[100,999]$ 中这些数的 $len(i)$ 都是固定的。于是我们对n 进行分段考虑。

接下来我们推导如何用矩阵转移这个 $dp$ 递推式

我们先对 $dp[i]$ 中的项进行转化，全部变成包含 $i-1$ 的式子。

dp[i]=10^{len(i)}\times dp[i-1]+(i-1)+1

则我们创建两个 $1\times 3$ 的向量 $I_{n},I_{n-1}$ 以及一个 $3\times 3$ 的矩阵 $A$ ，

则有以下关系：

I_{n+1}^T=A\times I_{n}^T

由上述状态转移方程有：

I_{n+1}^T= \begin{bmatrix} dp[n]\\ n\\ 1 \end{bmatrix}

以及：

I_{n}^T=\begin{bmatrix} dp[n-1]\\ n-1\\ 1 \end{bmatrix}

故有：

\begin{bmatrix} dp[n]\\ n\\ 1 \end{bmatrix}= \begin{bmatrix} a &b &c\\ d &e &f\\ g &h &i \end{bmatrix} \times \begin{bmatrix} dp[n-1]\\ n-1\\ 1 \end{bmatrix}

由于 $dp[i]=10^{len(i)}\times dp[i-1]+(i-1)+1$ ,所以有 $a=10^{len(i)},b=1,c=1$ ，

$n=(n-1)+1$ ，所以 $d=0,e=1,f=1$ ，

1=1，所以 $g=h=0,i=1$ 。

可以得出矩阵 $A$ ：

\begin{bmatrix} 10^{len(i)} &1 &1\\ 0 &1 &1\\ 0 &0 &1\\ \end{bmatrix}

由此可以得到：

\begin{bmatrix} dp[n]\\ n\\ 1 \end{bmatrix}= \begin{bmatrix} 10^{len(i)} &1 &1\\ 0 &1 &1\\ 0 &0 &1\\ \end{bmatrix} \times \begin{bmatrix} dp[n-1]\\ n-1\\ 1 \end{bmatrix}

最后再对 $len(i)$ 的长度进行 $dp$ 状态转移矩阵快速幂即可：

Code.:

struct Mat{
    i64 mat[4][4];
    Mat(){
        memset(mat,0,sizeof(mat));//矩阵初始化为0
    }
    
    void init_dp(i64 sz){//转移矩阵初始化
        mat[1][2]=1;
        mat[1][3]=1;
        mat[2][2]=1;
        mat[2][3]=1;
        mat[3][3]=1;
        mat[1][1]=sz;
        mat[1][1]%=m;
    }
    void init_b(){//初始化为单位矩阵
        mat[1][1]=1;
        mat[2][2]=1;
        mat[3][3]=1;
    }
    Mat operator * (const Mat &other){//结构体定义矩阵乘法
        Mat temp;
        for(int i=1;i<=3;i++){
            for(int j=1;j<=3;j++){
                i128 res=0;
                for(int k=1;k<=3;k++){
                    res+=mat[i][k]*other.mat[k][j];
                    res%=m;
                }
                temp.mat[i][j]=res;
            }
        }
        return temp;
    }
};
Mat binpow(Mat A,i64 x,i64 mod){//矩阵快速幂
    Mat I;
    I.init_b();
    while(x){
        if(x&1) I=I*A;
        A=A*A;
        x>>=1;
    }
    return I;
}
void solve(){
    cin>>n>>m;
    i64 f=0;
    for(i128 i=1,L=1;L<=n;i++,L*=10){
        i128 R=min((i128)L*10-1,(i128)n);
        i128 d=R-L+1;
        Mat base;
        base.init_dp(binpow(10,i,m));
        Mat temp=binpow(base,d,m);
        i64 nf=(temp.mat[1][1]*f+temp.mat[1][2]*((L-1)%m)+temp.mat[1][3])%m;
        f=nf;
    }

    cout<<f<<endl;
}

记得开__int128，不然会溢出超时。

END.

矩阵加速DP转移(自用）

矩阵加速DP转移

前置知识：矩阵快速幂

矩阵加速

例题

题目描述