矩阵快速幂优化DP

· · 算法·理论

可能更好的阅读体验

0. 前言

蒟蒻做题比较少，在做过的题中选出了一些经典的例题与技巧帮助大家，这篇文章也只是我个人的一个经验总结，希望能帮助到后人学习。

1.矩阵小芝士

矩阵优化是干啥的啊？

有的时候，你会发现你设计了一个极好的 DP 状态，没有后效性，没有重叠，你很开心，你去看数据范围就会炸掉！你死活都想不出来怎么优化，感觉去掉这个状态之后就感觉 “不完美” 了，让后点开题解，发现一堆密密麻麻的数学公式和矩阵，开心的点进去，郁闷的叉掉，那么怎么解决这种郁闷心情呢？当然就是学矩阵优化啦。

好端端的那么优美的递推式子我们为什么要用矩阵来转移呢？答案很简单，因为矩阵乘法有结合律，可以快速幂！当转移式子固定的时候，我们可以通过快速幂，设置好初始矩阵和转移矩阵，通过矩阵快速幂就能轻轻松松的将复杂度变成 \log，是不是很好？

我们介绍一下矩阵乘法：

设 A=(a_{i,j})_{m\times n},B=(b_{i,j})_{n\times s}，定义 C=A\times B=(c_{i,j})_{m\times s}，其中:

c_{i,j}=\sum_{k=1}^na_{i,k}b_{k,j}

即 A 的第 i 行与 B 的第 j 列相乘得到 C 的 (i,j) 元素。

向量可以视为 1\times n 或 n\times 1 的矩阵，类似定义向量和矩阵的乘法。

矩阵乘法满足结合律和分配律，但不满足交换律。

不只是普通乘法，我们还有广义矩阵乘法！

一些约定：

若 \otimes 满足交换律、结合律，\otimes 对 \oplus 有分配律，那么我们就有广义矩阵乘法：

(A\times B)_{i,j}=\bigoplus_{k}(A_{i,k}\otimes B_{j,k})

广义矩阵乘法也同样满足普通矩阵乘法的结合律、分配律。

例如 \oplus 是 \min 或 \max，\otimes 是 +。

广义上说，只要我们内层的运算具有结合律，外层的运算对于内层运算有分配律，并且状态转移递推恒定，并且发现转移来的状态数少但转移次数极大，我们就可以考虑矩阵优化 DP。

严谨的讲，其实矩阵优化的本质就是线性递推 DP 方程，也就是从 i 推到 i+1 的状态。比如说 f(i,j) 我们对 j 一维优化，通过转移矩阵我们可以将 f(i,j) \rightarrow f(i,j+1)。

设矩阵规模为 n\times n，转移阶段为 m，则复杂度为 O(n^3 \log m)，瓶颈在矩阵快速幂。

我们这里给出一个封装好的矩阵，其中 MN 需要额外定义常量。

struct Matrix{
    int mat[MN][MN];

    Matrix(int x=0){
        memset(mat,0,sizeof(mat));
        if(!x) return;
        for(int i=0;i<MN;i++) mat[i][i]=x;
    }

    Matrix operator*(const Matrix x)const{
        Matrix ret;
        for(int i=0;i<MN;i++){
            for(int j=0;j<MN;j++){
                for(int k=0;k<MN;k++){
                    ret.mat[i][j]+=mat[i][k]*x.mat[k][j];
                }
            }
        }
        return ret;
    }

};

Matrix ksm(Matrix a,int b){
    Matrix ret(1);
    while(b){
        if(b&1) ret=ret*a;
        a=a*a;
        b>>=1;
    }
    return ret;
}

2.矩阵优化操作手册

Luogu P1962 斐波那契数列

斐波那契数列是满足如下性质的一个数列：

f_n = \left\{\begin{aligned} 1 \space (n \le 2) \\ f_{n-1}+f_{n-2} \space (n\ge 3) \end{aligned}\right.

请你求出 f_n \bmod 10^9 + 7 的值。

这个题看起来还挺简单的，不就是斐波那契数列递推吗，直接 O(n) 乱搞就可以了。

但是 n 太大了怎么办？

我们发现，如果你想递推出 f_n 你需要用到 f_{n-1},f_{n-2} 的数据，而且这个也类似于一个 “转移”，而且不难发现每一次只用到上面两个的数据，转移来的状态极少，我们就可以考虑矩阵优化喽。

这里解释一下，转移来的状态指的是能够给当前状态贡献答案的状态。

通过简单的例题我们来介绍以下矩阵优化的 “操作手册”。

写出转移方程，判断优化类型

对于一般的矩阵优化方程，它的转移方程一般是不难写出来的，但是你会通常会发现其中一个维度数量极其大，转移十分困难，空间也十分困难，而且一个特征就是转移的式子很简单，近似于线性递推而且转移来的状态数量极少，这种就是矩阵优化的鲜明特征。

这里我们的方程就是：

f_n = \left\{\begin{aligned} 1 \space (n \le 2) \\ f_{n-1}+f_{n-2} \space (n\ge 3) \end{aligned}\right.

然而有的时候你可能看不出来维度数量极其大，这种我的建议就是你在设计 DP 的时候要写全，就是方程写出来，初始化设置好，计算好空间。

确定优化维度

我们要确定我们是要对哪一维度，普遍都是对那些数极其大的状态维度进行优化，有的时候可能不太一样，需要自行判断。

这里的优化维度显然就是 f_n 的 n。

根据转移方程需要的量，确定初始矩阵

确定初始矩阵实际上就是我们到底在转移需要什么数据才能转移到 i，这些数据可以通过转移方程来知道。

注意到，我们的转移方程需要的是前两个数，于是我们的初始矩阵可以这么设置：

f_n & f_n-1 \end{bmatrix}

你这不对吧，f_{n-2} 呢？

事实上，你应当回忆第一节后面我们讲的，矩阵优化的本质就是线性递推 DP 方程。我们从 f_n 推到 f_{n+1} 需要什么数据呢？回看转移方程，我们需要 f_{n} 和 f_{n-1}，所以我们在初始矩阵这样设置。

那么初始设置就是 [1,1]，与方程对应。

初始矩阵的实质？它是动态规划中递推起点的状态值所构成的矩阵，必须包含递推所需的全部初始状态，只有包含这些初始状态才能推到后面的状态。

一般来说，我们的初始矩阵大多是是一个向量的形式存在，我们的初始矩阵其实就是用来存状态的。

大多数时候初始矩阵的如果不填则默认为 0，但有的时候因为运算的限制我们需要填写负无穷（例如出现负数的时候或者取 \max 操作）

设计转移矩阵

难点来了，我们考虑如何设计转移矩阵。

我们回忆一下矩阵乘法的操作：

图画的不好，请见谅，但大体上是这个操作，我们考虑我们的转移方程该怎么操作。

首先明确目标，我们是要从 f_n 递推到 f_{n+1}。我们目前有转移到 f_{n+1} 的必须数据——前两项。

接下来我们确定转移后的数据，转移后我们要递推 f_{n+2}，就需要 f_{n+1},f_n 的数据，而原来矩阵的 f_{n-1} 的数据以及没用了，我们就可以舍弃，那么我们列出初始矩阵和转移后的目标矩阵。

接下来是填空时间！请读者自行填空。若感到吃力可以看上面我们矩阵运算的示意图。答案下面会给出

我们总结一下这个步骤：

  1. 确定转移矩阵大小。
  2. 确定初始矩阵转移后的目标矩阵。
  3. 根据转移方程和目标矩阵，列出转移矩阵。
  4. 利用转移矩阵模拟运算后是否能得出目标矩阵的结果。

答案如下：

确定转移顺序，确定快速幂的幂数

我们不妨设初始矩阵为 A，转移矩阵为 G。

那么求 f_n 的过程就是 A\times G \times G \times \dots \times G。

注意到一开始说矩阵乘法是有结合律的，考虑直接把后面结合起来，注意到是一个幂的形式，于是我们就可以快速幂了！

那为什么需要这一步呢，是因为有的时候转移过程中因为题目条件矩阵转移顺序有所变化，所以我们可能需要中断转移改变转移矩阵，这里下面的例题会讲到。

在这里我们需要确定矩阵快速幂的幂数，注意到我们初始矩阵是 [f_2,f_1] （不能用 f_0 因为就没有第 0 项），那么我们总共需要递推 n-2 次，所以快速幂的幂数就是 n-2。

根据初始矩阵设计，确定答案统计范围

这个一般需要你通过初始矩阵设计来确定，可能是矩阵的某一项，或者矩阵一个范围。

写代码，注意细节。

写代码的时候尤其要注意矩阵初始化的赋值，到底是 0 还是初始化为正无穷还是负无穷。注意矩阵的数据大小选择合适的存储类型。

#include<bits/stdc++.h>
#define ll long long
using namespace std;
const int MOD=1e9+7,MT=5;
struct matrix{
    ll mt[MT][MT]{};
    matrix operator *(const matrix &x){
        matrix ret;
        for(int i=1;i<=2;i++){
            for(int j=1;j<=2;j++){
                for(int k=1;k<=2;k++){
                    ret.mt[i][j]+=mt[i][k]*x.mt[k][j]%MOD;
                    // ret.mt[i][j]%=MOD;
                }
            }
        }
        return ret;
    }
};
ll n;

matrix qp(matrix x,matrix ans,ll k){
    while (k)
    {
        if(k&1){
            ans=ans*x;
        }
        x=x*x;
        k>>=1;
    }
    return ans;
}
matrix ans,base;
int main(){
    cin>>n;
    if(n<2){
        cout<<n;
        return 0;
    }
    base.mt[1][1]=base.mt[1][2]=base.mt[2][1]=1;
    ans.mt[1][1]=ans.mt[1][2]=1;
    cout<<qp(base,ans,n-2).mt[1][1]%MOD;
    return 0;
}

3. 例题与技巧

接下来是例题时间，对于每一个例题我们后面都会有对应的技巧进行讲解。

USACO07NOV Cow Relays G——Floyd倍增转移

给定一张 T 条边的无向连通图，求从 S 到 E 经过 N 条边的最短路长度。

你真的理解 Floyd 了吗？

我们回忆图的邻接矩阵，邻接矩阵的本质是矩阵，它原始表示图两点之间经过 1 条边的路径数量。

我们考虑 Floyd 算法本质上是在干什么，对于 3 层循环 (k,i,j)，其中 k 我们在枚举中间点，我们来看 Floyd 的转移：mp[i][j]=min(mp[i][j],mp[i][k]+mp[k][j]);。你会发现，这有点类似于广义矩阵乘法啊！

事实上，Floyd 枚举中间点的转移，事实上对于原来的邻接矩阵来说，从经过 1 条边到 2 条边，从经过 2 条边到经过 3 条边，如此下去将所有边遍历完就能求出多源最短路了！

以下我们改写以下方程：

c[i][j]=\min(c[i][j],a[i][k]+b[k][j])

根据 Floyd，a 矩阵是经过 x 条边的最短路，而 b 为经过 y 条边的最短路，那么 c 矩阵就是经过 x+y 的最短路，那么我们根据初始输入的数组，一开始是经过 1 条边的，那么转移 n-1 次就是我们想要的答案了，我们显然不能暴力转移，注意到满足广义矩阵乘法进行快速幂即可，注意实现细节。

#include<bits/stdc++.h>
using namespace std;
constexpr int MN=520;
int n,t,s,e,tot;
struct Matrix{
    int mt[MN][MN];
    Matrix operator*(const Matrix &x)const{
        Matrix ans;
        memset(ans.mt,0x3f,sizeof(ans.mt));
        for(int k=1;k<=tot;k++){
            for(int i=1;i<=tot;i++){
                for(int j=1;j<=tot;j++){
                    ans.mt[i][j]=min(ans.mt[i][j],mt[i][k]+x.mt[k][j]);
                }
            }
        }
        return ans;
    }
}dis,ans;
unordered_map<int,int> ump;

Matrix ksm(){
    n--;
    ans=dis;
    while(n){
        if(n&1){
            ans=ans*dis;
        }
        dis=dis*dis;
        n>>=1;
    }
    return ans;
}

int main(){
    cin>>n>>t>>s>>e;
    memset(dis.mt,0x3f,sizeof(dis.mt));
    for(int i=1;i<=t;i++){
        int u,v,w;
        cin>>w>>u>>v;
        if(!ump[u]) ump[u]=++tot;
        if(!ump[v]) ump[v]=++tot;
        dis.mt[ump[u]][ump[v]]=dis.mt[ump[v]][ump[u]]=w;
    }
    cout<<ksm().mt[ump[s]][ump[e]];

    return 0;
}

经典图例题

给定一个 n 个点，m 条边的有向图，对于两个点之间要么不连边，要么边权为 1，求从起点出发长度为 k 到达终点的路径方案数（n \le 500,m\le 250，k 极大）。

先不考虑数据范围，我们可以用 DP 来解决。

设 f(i,j) 表示到第 i 个点，路径长度为 j 的路径条数，注意到边权为 1 我们可以考虑直接通过相邻边 BFS 转移即可，时间复杂度 O(mk)。

注意到 k 极大但 m 很小，也就是转移过来的状态数量极少，自然想到矩阵优化，那么我们根据什么列转移矩阵呢？我们可以考虑邻接矩阵，我们利用邻接矩阵类似 Floyd 的方法来进行转移，时间复杂度 O(n^3 \log k)。

SCOI2009迷路——边权拆点

给定一个 n 个点，m 条边的有向图，边权为 w，求从 1 号点出发长度为 k 到达 n 号点的路径方案数（n \le 10,w \le 9，k\le 10^9）。

注意到这里边权为 w。我们矩阵的一个要求就是线性递推，从 k-1 \rightarrow k，但是这里因为有边权就不满足了，如何处理？

这里我们用到一个技巧叫做边权拆点，矩阵的要求就是线性递推，那么为了满足线性递推的要求，我们把一条边拆成一个一个边权为 1 的小边，同时我们把一个节点拆成 9 个点，分别表示周围边权为 1 到周围边权为 9。对于拆点内部，我们让 (u,i)\rightarrow (u,i+1) 连边；对于拆点外部，我们让 (u,w) \rightarrow (v,1)。 例如 u\rightarrow v,w=4 这条边，我们就是将 (u,4) 连向 (v,1) 号节点。如此，这样构建能够保证当前图和原图等价，于是直接跑矩阵加速即可。

使用这个技巧的的前提是边权很小，这样我们才能用边权拆点。

拆点实际上就是将原图换为一个于此等价的图，使得题目中的特殊性质被满足，从而达到简化题目的目的。

我们其实也可以拆边，把边拆为点，但是我们应当在尽量满足题意的情况下尽量减少点的数量，从而减低复杂度，时间复杂度为 O((9n)^3 \log k)。

#include<bits/stdc++.h>
using namespace std;
constexpr int MN=100,MOD=2009;
int n,t,tot,idx[MN][10];

struct Matrix{
    int mat[MN][MN];

    Matrix(int x=0){
        memset(mat,0,sizeof(mat));
        if(!x) return;
        for(int i=0;i<MN;i++) mat[i][i]=x; 
    }

    Matrix operator*(const Matrix &x)const{
        Matrix ret;
        for(int i=0;i<MN;i++){
            for(int j=0;j<MN;j++){
                for(int k=0;k<MN;k++){
                    ret.mat[i][j]+=mat[i][k]*x.mat[k][j];
                    ret.mat[i][j]%=MOD;
                }
            }
        }
        return ret;
    }

}A,G;

Matrix ksm(Matrix a,int b){
    Matrix ret(1);
    while(b){
        if(b&1) ret=ret*a;
        a=a*a;
        b>>=1;
    }
    return ret;
}

void initG(){
    for(int i=1;i<=n;i++){
        for(int j=1;j<=9;j++){
            idx[i][j]=++tot;
        }
    }
    for(int i=1;i<=n;i++){
        for(int j=1;j<9;j++){
            G.mat[idx[i][j]][idx[i][j+1]]=1;
        }
    }
}

int main(){
    cin>>n>>t;
    initG();
    for(int i=1;i<=n;i++){
        for(int j=1;j<=n;j++){
            char c;
            int num;
            cin>>c;
            num=c-'0';
            if(num>0){
                G.mat[idx[i][num]][idx[j][1]]=1;
            }
        }
    }
    A=ksm(G,t);
    cout<<A.mat[1][idx[n][1]]%2009;
    return 0;
}

SDOI2009 HH去散步——点边转化，超级源点

给定 n 个点，m 条无向边的图，边权均为 1，求从 S\rightarrow T 的路径上有多少长度为 t 的路径，满足上一步走过的路下一步不能重复走。答案对 45989 取模。

考虑DP，设 f(i,j) 表示走到第 i 个点，路径长度为 j 的路径条数，显然有转移方程：

\begin{aligned} f(u,j)=\sum\limits_{v\in son(u)} f(v,j-1) \end{aligned}

那你这方程也不对啊，题目说了满足上一步走过的路下一步不能重复走。那么接下来怎么在图上转移？

卡住了……但是我们发现边数极小！当然我们需要构造一个新图使得能够满足这个限制，考虑怎么构造？

注意到边数极小，我们可以利用点边互换的技巧！点边互换的核心思想是：把原图中的某些点看作边，或者把原图中的某些边看作点 ，从而构造一个新的图来满足题目中的限制。

边可以转为点，而点也可以转为边，这样的代价就是我们必须增加点的数量或边的数量。

那么我们修改一下转移方程，转移方程基本不变，但是转移有区别。

f(u,j)=\sum\limits_{v\in to(u)} f(v,j-1)

其中 to(u) 表示能够到达 u 点的边集，那么这个边集怎么搞呢，我们可以考虑割点！把每个点拆成相邻边数个点，第 i 个点代表第 i 条相邻的边走过来得点，让后该点向除第 i 条边以外所连的点连一条单向边即可构造出新图。

但是 S 节点可能向周围任意一个点前进，不好在矩阵掌控，不妨考虑超级源点 0 号点，0 号点向 S 连边，只需要把 f_0 初始化为 1 即可跑矩阵优化。

对于每一个点，期望连边 2 条，都会割出 2 个点，最多 120 个点，时间复杂度 (m^3 \log k)。

割点其实就是将每个点的不同状态放大化，使其“人格分裂”，从而将【节点的转移】转化为【状态的转移】，从而达到实现满足特殊条件约束的目的。 ——Ydtz 的教程

HAOI2015 数字串拆分——指数高精？发掘性质！

自己看题吧。

注意到 f 可以 DP 求出，根据题意列出转移方程。

f(i)=\sum\limits_{j=max(i-m,0)}^{i-1} f(j)

初始 f(0)=1。

注意到 m 极小，而 |s_0| 及其大，考虑矩阵优化，注意到这个和斐波那契数列的转移有点类似，不过要的是前 m 个数。

但是你会发现这个 g 的划分及其奇怪，并且转移的 10^{500} 太大了要指数高精，而且暴力划分 g 最大会有 2^{|s_0|} 种。

考虑分析 g 的划分形式，首先矩阵满足一个性质：

A^{n+m}=A^n \times A^m

正确性是很显然的，例如 g(123) 可以分解为：

\begin{aligned} g(123) & =f(1+2+3)+f(12+3)+f(1+23)+f(123) \\ & = f(1+2) \times 3+f(12) \times f(3) +f(1) \times f(23)+f(123) \end{aligned}

考虑按数位递推求解 g，设 g_i' 表示对于 s 中前 i 位数字的答案 g，考虑算 g(123) 时 g(12),g(1) 已经算过，考虑改写有：

\begin{aligned} g(123) & = f(1+2) \times 3+f(12) \times f(3) +f(1) \times f(23)+f(123) \\ & = f(123) + g'_{1} \times f(23) +g'_{2} \times f(3) \end{aligned}

答案即为 g'_n。

问题转化为如何求出一个后缀字符串的 f 值，考虑 f'_{i,j} 表示 s 第 i 位数字第 j 位数字组成的 f 的转移矩阵，通过 f'_{i,j}=f(10^{n-i})^{s_{i}}\times f'_{i+1,j} 即可转移，预处理 f(10^i) 即可，其中 s_i 表示代表 s 中第 i 个数字的矩阵。

当我们发现指数需要高精，心中默念 “这题绝对不是让我指数高精，一定题目中有什么其他的特性在里面，如果山穷水尽可能我们还有特殊快速幂” 关于这个特殊快速幂我们下文会提到。

#include<bits/stdc++.h>
#define int long long
using namespace std;
constexpr int MN=5,MOD=998244353;
int m,len,ans;
string s;

struct Matrix{
    int mat[MN][MN];

    Matrix(int x=0){
        memset(mat,0,sizeof(mat));
        if(x==0) return;
        for(int i=0;i<MN;i++){
            mat[i][i]=x;
        }
    }

    void init(int x){
        memset(mat,0,sizeof(mat));
        if(x==0) return;
        for(int i=0;i<MN;i++){
            mat[i][i]=x;
        }
    }

    void initf(){
        for(int i=0;i<m;i++) mat[i][m-1]=1;
        for(int i=1;i<m;i++) mat[i][i-1]=1;
    }

    Matrix operator*(const Matrix &x)const{
        Matrix ret;
        for(int i=0;i<MN;i++){
            for(int j=0;j<MN;j++){
                for(int k=0;k<MN;k++){
                    ret.mat[i][j]+=mat[i][k]*x.mat[k][j]%MOD;
                    ret.mat[i][j]%=MOD;
                }
            }
        }
        return ret;
    }

    Matrix operator+(const Matrix &x)const{
        Matrix ret;
         for(int i=0;i<MN;i++){
            for(int j=0;j<MN;j++){
                ret.mat[i][j]=(mat[i][j]+x.mat[i][j])%MOD;
            }
        }
        return ret;
    }

}pw[520],f[520][520],g[520];

Matrix ksm(Matrix a,int b){
    Matrix ret(1);
    while(b>0){
        if(b&1) ret=ret*a;
        a=a*a;
        b>>=1;
    }
    return ret;
}

signed main(){
    cin>>s>>m;
    len=s.length();
    s=' '+s;
    pw[0].initf();
    for(int i=1;i<len;i++) pw[i]=ksm(pw[i-1],10);
    for(int i=1;i<=len;i++){
        for(int j=i;j>=1;j--){
            if(i!=j){
                f[j][i]=f[j+1][i]*ksm(pw[i-j],s[j]-'0');
            }else{
                f[j][i]=ksm(pw[0],s[j]-'0');
            }
        }
    }
    g[0].init(1);
    for(int i=1;i<=len;i++){
        for(int j=0;j<i;j++){
            g[i]=g[i]+(g[j]*f[j+1][i]);
        }
    }
    for(int i=0;i<m;i++) ans=(ans+g[len].mat[0][i])%MOD;
    cout<<ans;
    return 0;
}

NOI2013 矩阵游戏——指数高精与费马小定理，10进制快速幂。

F[1, 1] &= 1 \\ F[i, j] &=a\times F[i, j-1]+b, &j\neq 1 \\ F[i, 1] &=c\times F[i-1, m]+d, &i\neq 1 \\ \end{aligned}

不难发现是线性递推，考虑设置两个转移方程 M1,M2 分别表示转移列的和转移行的。

转移方程已经写出来了，我们考虑怎么设置初始转移，我们不难发现这个矩阵转移常数项怎么处理，对于常数项的处理我们一般是在初始矩阵中加一位列单独设置为 1，让后在之后的转移一直设置为 1，如果需要常数项我们就将对应位置改成 b 或 d 即可。

矩阵如下,列转移：

f & 1 \end{bmatrix} \times \begin{bmatrix} a & 0\\ b & 1 \end{bmatrix} = \begin{bmatrix} f' & 1 \end{bmatrix}

行转移：

f & 1 \end{bmatrix} \times \begin{bmatrix} c & 0\\ d & 1 \end{bmatrix} = \begin{bmatrix} f' & 1 \end{bmatrix}

对于转移实际上就是 A \times ({M1}^{m-1}\times M2)^{n-1} \times {M1}^{m-1}。

注意到 n,m 的范围及其离谱，\log_2 都搞不了。

等会 \log_2 不行其他的 \log_x 是不是可以，我们可以考虑 10进制快速幂，让后就可以了！

显然我没写 10进制快速幂，心中默念是不是有没有什么性质。

注意到这个矩阵对于列的转移，其中 a 总共被乘上了 a^{m-1} 次，同理与 c^{n-1}，而注意到模数为质数，可以考虑费马小定理取模，模上 \varphi(MOD) 即可。

注意到，当 a=c=1 时不能用费马小定理，考虑直接矩阵快速幂的幂取模原来模数。

#include<bits/stdc++.h>
#define int __int128
using namespace std;
constexpr int MN=15,MOD=1e9+7;
int n,m,a,b,c,d,base0,base1;
string sn,sm;

struct Matrix{
    int mat[MN][MN];

    // INIT THE MATRIX
    Matrix(int x=0){
        memset(mat,0,sizeof(mat));
        for(int i=0;i<MN;i++) mat[i][i]=x;
    }

    Matrix operator*(const Matrix &x)const{
        Matrix ret;
        for(int i=0;i<MN;i++){
            for(int j=0;j<MN;j++){
                for(int k=0;k<MN;k++){
                    ret.mat[i][j]+=mat[i][k]*x.mat[k][j];
                    ret.mat[i][j]%=MOD;
                }
            }
        }
        return ret;
    }

}M1,M2;

Matrix ksm(Matrix A,int B){

    Matrix ret(1);
    while(B>0){
        if(B&1) ret=ret*A;
        A=A*A;
        B>>=1;
    }
    return ret;
}

namespace ly
{
    namespace IO
    {
        #ifndef LOCAL
            constexpr auto maxn=1<<20;
            char in[maxn],out[maxn],*p1=in,*p2=in,*p3=out;
            #define getchar() (p1==p2&&(p2=(p1=in)+fread(in,1,maxn,stdin),p1==p2)?EOF:*p1++)
            #define flush() (fwrite(out,1,p3-out,stdout))
            #define putchar(x) (p3==out+maxn&&(flush(),p3=out),*p3++=(x))
            class Flush{public:~Flush(){flush();}}_;
        #endif
        namespace usr
        {
            template<typename type>
            inline type read(type &x)
            {
                x=0;bool flag(0);char ch=getchar();
                while(!isdigit(ch)) flag^=ch=='-',ch=getchar();
                while(isdigit(ch)) x=(x<<1)+(x<<3)+(ch^48),ch=getchar();
                return flag?x=-x:x;
            }
            template<typename type>
            inline void write(type x)
            {
                x<0?x=-x,putchar('-'):0;
                static short Stack[50],top(0);
                do Stack[++top]=x%10,x/=10;while(x);
                while(top) putchar(Stack[top--]|48);
            }
            inline char read(char &x){do x=getchar();while(isspace(x));return x;}
            inline char write(const char &x){return putchar(x);}
            inline void read(char *x){static char ch;read(ch);do *(x++)=ch;while(!isspace(ch=getchar())&&~ch);}
            template<typename type>inline void write(type *x){while(*x)putchar(*(x++));}
            inline void read(string &x){static char ch;read(ch),x.clear();do x+=ch;while(!isspace(ch=getchar())&&~ch);}
            inline void write(const string &x){for(int i=0,len=x.length();i<len;++i)putchar(x[i]);}
            template<typename type,typename...T>inline void read(type &x,T&...y){read(x),read(y...);}
            template<typename type,typename...T>
            inline void write(const type &x,const T&...y){write(x),putchar(' '),write(y...),sizeof...(y)^1?0:putchar('\n');}
            template<typename type>
            inline void put(const type &x,bool flag=1){write(x),flag?putchar('\n'):putchar(' ');}
        }
        #ifndef LOCAL
            #undef getchar
            #undef flush
            #undef putchar
        #endif
    }using namespace IO::usr;
}using namespace ly::IO::usr;

signed main(){
    read(sn,sm,a,b,c,d);
    base0=(a==1?MOD:MOD-1);
    for(auto p:sm){
        m=((m<<3)+(m<<1)+(p^48))%(base0);
    }
    Matrix ans;
    M1.mat[0][0]=a,M1.mat[0][1]=b,M1.mat[1][1]=1;
    M2.mat[0][0]=c,M2.mat[0][1]=d,M2.mat[1][1]=1;
    ans.mat[0][0]=ans.mat[1][0]=1;
    Matrix d=ksm(M1,(m+base0-1)%base0)*M2;
    int base1;
    if(d.mat[0][0]==1) base1=MOD;
    else base1=MOD-1;
    for(auto p:sn){
        n=((n<<3)+(n<<1)+(p^48))%(base1);
    }
    ans=ksm(d,(n+base1-1)%base1)*ksm(M1,(m+base0-1)%base0)*ans;
    put(ans.mat[0][0]);
    return 0;
}

CF576D Flights for Regular Customers——强制中断，bitset优化

又是特殊限制，我们还是设 DP。

设 f(i,j) 表示在第 i 个点，在走过的边数为 j 的情况下是否能够到达（取值为 0 或 1），由 j-1 可以转移过来，并且矩阵味很重，转移是或的关系，可以考虑矩阵优化。

考虑无解的情况怎么做，不妨假设 1 号节点边都可以走，如果都可以走的情况下还是到不了那就 GG。

我们根据操作手册，发现在第五步就炸了，因为每一次 d_i 的更新都需要重新设置转移矩阵，考虑根据 d_i 的变化量进行快速幂，每一次中断跑多源 BFS 更新答案，让后就做完了。

我们不难发现 f 的取值只有 0 或 1，可以利用 bitset 优化，写的时候如下：

struct Matrix{
    bitset<MN> mat[MN];

    Matrix(int x=0){
        for(int i=0;i<MN;i++){
            for(int j=0;j<MN;j++){
                mat[i][j]=0;
            }
        }
        if(!x) return;
        for(int i=0;i<MN;i++) mat[i][i]=x;
    }

    Matrix operator*(const Matrix &x)const{
        Matrix ret;
        for(int i=0;i<MN;i++){
            for(int k=0;k<MN;k++){
                if(mat[i][k]){// 把j省去了
                    ret.mat[i]|=x.mat[k];
                }
            }
        }
        return ret;
    }

};

代码如下：

#include<bits/stdc++.h>
#define int long long
using namespace std;
constexpr int MN=250,INF=1e18;
struct Edge{
    int u,v,d;
}e[MN];
int n,m,ans,dis[MN];

struct Matrix{
    bitset<MN> mat[MN];

    Matrix(int x=0){
        for(int i=0;i<MN;i++){
            for(int j=0;j<MN;j++){
                mat[i][j]=0;
            }
        }
        if(!x) return;
        for(int i=0;i<MN;i++) mat[i][i]=x;
    }

    Matrix operator*(const Matrix &x)const{
        Matrix ret;
        for(int i=0;i<MN;i++){
            for(int k=0;k<MN;k++){
                if(mat[i][k]){
                    ret.mat[i]|=x.mat[k];
                }
            }
        }
        return ret;
    }

}a,G;

Matrix ksm(Matrix a,int b){
    Matrix ret(1);
    while(b>0){
        if(b&1) ret=ret*a;
        a=a*a;
        b>>=1;
    }
    return ret;
}

bool cmp(Edge x,Edge y){
    return x.d<y.d;
}

void bfs(){
    for(int i=1;i<=n;i++) dis[i]=INF;
    queue<int> q;
    for(int i=1;i<=n;i++){
        if(a.mat[1][i]) q.push(i),dis[i]=0;
    }
    while(!q.empty()){
        int f=q.front();
        cerr<<f<<'\n';
        q.pop();
        for(int i=1;i<=n;i++){
            if(G.mat[f][i]&&dis[i]==INF){
                dis[i]=dis[f]+1;
                q.push(i);
            }
        }
    }
}

signed main(){
    cin>>n>>m;
    for(int i=1;i<=m;i++){
        cin>>e[i].u>>e[i].v>>e[i].d;
    }
    sort(e+1,e+1+m,cmp);
    a.mat[1][1]=1,ans=INF,dis[n]=INF;
    for(int i=1,lst=0;i<=m;i++){
        if(ans<e[i].d) break;
        int dt=e[i].d-lst;
        lst=e[i].d;
        a=a*ksm(G,dt);
        G.mat[e[i].u][e[i].v]=1;
        if(i==m||e[i+1].d!=e[i].d) bfs();
        ans=min(ans,lst+dis[n]);
    }
    if(ans==INF) cout<<"Impossible";
    else cout<<ans;
    return 0;
}

NOI2020 美食家——强制中断，二进制分组快速幂！

给定一张 n 个点 m 条边的有向图，有重边，边有时间 w_i。每一个点有价值 c，走到点 i 可以获得 c_i 的价值。 初始时间为 0，你需要从起点 1 开始走，恰好在 T 时间回到起点 1。最终得到的价值为所有经过的点的价值 c 的和，点可以被重复经过而且价值也可以重复累加。 不能在任何一个点停留，也就是说 t 时间到达一个点 t+1 时间必须往其他点走。 现在有 k 个特殊时间点，形式为三个参数 (t_i,x_i,y_i)。表示恰好在 t_i 时间点到达点 x_i 时可以得到 y_i 的额外价值。 求最终能获得的最大价值和，若无法在 T 时间回到起点，输出 -1。

1\le n\le 50,n\le m \le 500,0 \le k \le 200 1\le t_{i}\le T \le 10^9,1\le w_{i}\le 5,1\le c_{i}\le 52501,1\le y_{i}\le 10^9

数据保证所有 t_i 互不相同，图一定联通。

一眼 DP（因为我也想不出来如何用贪心做），考虑设 f(i,j) 表示到达第 i 个点，当前时间为 j 的最大价值和，因为没有停留甚至都不用设置是否停留的状态太好啦。转移方程是显然的：

\begin{aligned} f(u,j)=\max\limits_{v=fa(u)}f(v,j-1)+w(u,v) \end{aligned}

这是转移方程，我们可以写成递推的形式，就不写出了。考虑特殊时间，直接用 map 存下来特判即可，这样我们就能拿到高贵的 40 分（环特判时间即可），对于性质 A 是一个大环直接瞎做即可，就有 50 分。

首先特别重要的一点，1\le w_{i} \le 5,n\le m \le 500。如果这个图是一个完全图则对应的边数也太少了，每个点期望连边的边数很小，也就是说 DP 方程转移过来的状态数量极小，并且 w_i 极小，有因为取 \max 和内部操作满足广义矩阵乘法，考虑矩阵优化 + 边权拆点。

不对啊，还是卡在了第五步，因为有特殊时间点的存在，我们可以想上面题一样强制中断快速幂，更新矩阵后再次快速幂，ok 啊，写完，交上去发现只有 75 分（拼好分版本）。

怎么回事呢，原来是 k 太炸裂了，这样的时间复杂度为 O(k\times (5n)^3 \log T)，只能过 k\le 10。

我们考虑怎么优化，发现每次中断转移后重新转移的 \Delta t，\Delta t 中有一些幂我们是重复在乘上的。

我们考虑对 \Delta t 进行二进制分解，首先我们写预处理转移矩阵 G 的幂 pw[i]，(G^1,G^2,G^4,\dots,G^{2^{29}})，预处理的时间复杂度就是 O((5n)^3 \log T)。让后我们每一次快速幂转移的时候，我们将 \Delta t 做二进制分解，将对应二进制数上的幂对应预处理的转移矩阵幂。让后我们让初始矩阵乘上这些 pw[i]，就可以了。

分析复杂度，不难发现我们的初始矩阵是一个 1 行 5n 列的向量，复杂度为 O(k \times (5n)^2 \log T)，可以通过。

类似于这种凑 k 的多次询问或者需要多次用到凑 k，如果求解二进制整数幂答案的时间复杂度较低，可以考虑倍增+二进制拆分求解。

注意转矩阵乘法后千万不要思维定势，因为这里是取 \max 操作，所有无用项和初始值为负无穷而非 0。

#include<bits/stdc++.h>
#define int long long
using namespace std;
constexpr int MN=260,NINF=-1e18;
struct Festival{
    int t,u,y;
}fst[MN];
int logt,n,m,T,K,tot,idx[MN][MN],c[MN];

struct Matrix{
    int mat[MN][MN];

    Matrix(int x=NINF){
        for(int i=0;i<MN;i++){
            for(int j=0;j<MN;j++){
                if(i==j) mat[i][j]=x;
                else mat[i][j]=NINF;
            }
        }
    }

    Matrix operator*(const Matrix &x)const{
        Matrix ret;
        for(int i=0;i<MN;i++){
            for(int j=0;j<MN;j++){
                for(int k=0;k<MN;k++){
                    if(mat[i][k]==NINF||x.mat[k][j]==NINF) continue;
                    ret.mat[i][j]=max(ret.mat[i][j],mat[i][k]+x.mat[k][j]);
                }
            }
        }
        return ret;
    }

}pw[55],A,B;

void initA(){
    for(int i=1;i<=n;i++){
        for(int j=1;j<=5;j++){
            idx[i][j]=++tot;
        }
    }
    for(int i=1;i<=n;i++){
        for(int j=1;j<5;j++){
            A.mat[idx[i][j]][idx[i][j+1]]=0;   
        }
    }
}

void initpw(Matrix x){
    pw[0]=x;
    for(int i=1;i<=logt;i++) pw[i]=pw[i-1]*pw[i-1];
}

bool cmp(Festival x,Festival y){
    return x.t<y.t;
}

Matrix ksmpw(Matrix a,int b){
    for(int i=0;i<=logt;i++){
        if((b>>i)&1){
            Matrix ret;
            for(int j=0;j<MN;j++){
                for(int k=0;k<MN;k++){
                    if(a.mat[1][k]==NINF||pw[i].mat[k][j]==NINF) continue;;
                    ret.mat[1][j]=max(ret.mat[1][j],a.mat[1][k]+pw[i].mat[k][j]);
                }
            }
            a=ret;
        }
    }
    return a;
}

signed main(){
    cin>>n>>m>>T>>K;
    logt=__lg(T);
    for(int i=1;i<=n;i++){
        cin>>c[i];
    }
    initA();
    for(int i=1;i<=m;i++){
        int u,v,w;
        cin>>u>>v>>w;
        A.mat[idx[u][w]][idx[v][1]]=c[v];
    }
    for(int i=1;i<=K;i++){
        cin>>fst[i].t>>fst[i].u>>fst[i].y;
    }
    sort(fst+1,fst+1+K,cmp);
    initpw(A);
    B.mat[1][idx[1][1]]=c[1];
    int lst=0;
    for(int i=1;i<=K;i++){
        B=ksmpw(B,fst[i].t-lst);
        if(B.mat[1][idx[fst[i].u][1]]!=NINF) B.mat[1][idx[fst[i].u][1]]+=fst[i].y;
        lst=fst[i].t;
    }
    if(lst!=T) B=ksmpw(B,T-lst);
    if(B.mat[1][idx[1][1]]==NINF) cout<<-1;
    else cout<<B.mat[1][idx[1][1]];
    return 0;
}

POI 2015 WYC——倍增二分思想

给定一个 n 个点，m 条边的有向图，无自环有重边，边权为 w，求第 k 小路径长度（1\le n \le 40,1\le m\le 1000,1 \le w \le 3，k\le 10^{18}）。

A* ? 痴人说梦啊。

注意到又是经典的 n,w 极小，k 极大，稍微分析以下满足矩阵的特点，考虑矩阵递推 + 边权拆点，因为有重边所以不能单纯的矩阵设为 1，应当为G.mat[idx[u][c]][idx[v][1]]++;。但是我们怎么统计一个状态有多少路径呢？我们发现可以利用超级源点的思想，0 号点向他们连边，边权为 0，这样就能够统计统计 0 号点的路径得到全局答案了。

注意到我们需要找第 k 小，也就是凑 k，我们就可以利用前面提到的思想。注意到这个是一个倍增的形式，我们可以利用倍增二分的技巧找 k，时间复杂度就是 O(n^3 \log k)。

注意开 __int128:

#include<bits/stdc++.h>
#define int __int128
using namespace std;
constexpr int MN=150;
int n,m,k,lgk,tot,idx[MN][5];

struct Matrix{
    int mat[MN][MN];

    Matrix(int x=0){
        memset(mat,0,sizeof(mat));
        if(!x) return;
        for(int i=0;i<MN;i++) mat[i][i]=x;
    }

    Matrix operator*(const Matrix x)const{
        Matrix ret;
        for(int i=0;i<MN;i++){
            for(int j=0;j<MN;j++){
                for(int k=0;k<MN;k++){
                    ret.mat[i][j]+=mat[i][k]*x.mat[k][j];
                }
            }
        }
        return ret;
    }

}a,G,pw[70];

namespace ly
{
    namespace IO
    {
        #ifndef LOCAL
            constexpr auto maxn=1<<20;
            char in[maxn],out[maxn],*p1=in,*p2=in,*p3=out;
            #define getchar() (p1==p2&&(p2=(p1=in)+fread(in,1,maxn,stdin),p1==p2)?EOF:*p1++)
            #define flush() (fwrite(out,1,p3-out,stdout))
            #define putchar(x) (p3==out+maxn&&(flush(),p3=out),*p3++=(x))
            class Flush{public:~Flush(){flush();}}_;
        #endif
        namespace usr
        {
            template<typename type>
            inline type read(type &x)
            {
                x=0;bool flag(0);char ch=getchar();
                while(!isdigit(ch)) flag^=ch=='-',ch=getchar();
                while(isdigit(ch)) x=(x<<1)+(x<<3)+(ch^48),ch=getchar();
                return flag?x=-x:x;
            }
            template<typename type>
            inline void write(type x)
            {
                x<0?x=-x,putchar('-'):0;
                static short Stack[50],top(0);
                do Stack[++top]=x%10,x/=10;while(x);
                while(top) putchar(Stack[top--]|48);
            }
            inline char read(char &x){do x=getchar();while(isspace(x));return x;}
            inline char write(const char &x){return putchar(x);}
            inline void read(char *x){static char ch;read(ch);do *(x++)=ch;while(!isspace(ch=getchar())&&~ch);}
            template<typename type>inline void write(type *x){while(*x)putchar(*(x++));}
            inline void read(string &x){static char ch;read(ch),x.clear();do x+=ch;while(!isspace(ch=getchar())&&~ch);}
            inline void write(const string &x){for(int i=0,len=x.length();i<len;++i)putchar(x[i]);}
            template<typename type,typename...T>inline void read(type &x,T&...y){read(x),read(y...);}
            template<typename type,typename...T>
            inline void write(const type &x,const T&...y){write(x),putchar(' '),write(y...),sizeof...(y)^1?0:putchar('\n');}
            template<typename type>
            inline void put(const type &x,bool flag=1){write(x),flag?putchar('\n'):putchar(' ');}
        }
        #ifndef LOCAL
            #undef getchar
            #undef flush
            #undef putchar
        #endif
    }using namespace IO::usr;
}using namespace ly::IO::usr;

void initG(){
    for(int i=1;i<=n;i++){
        for(int j=1;j<=3;j++) idx[i][j]=++tot;
    }
    for(int i=1;i<=n;i++){
        for(int j=1;j<3;j++){
            G.mat[idx[i][j]][idx[i][j+1]]++;

        }
        G.mat[idx[i][1]][0]++;
        a.mat[0][idx[i][1]]++;
    }
    G.mat[0][0]++;
}

signed main(){
    read(n,m,k);
    initG();
    for(int i=1;i<=m;i++){
        int u,v,c;
        read(u,v,c);
        G.mat[idx[u][c]][idx[v][1]]++;
    }
    int i;
    pw[0]=G;
    for(i=1;;i++){
        pw[i]=pw[i-1]*pw[i-1];
        Matrix ret=a*pw[i];
        if(ret.mat[0][0]-n>=k) break;
        if(i>62){
            cout<<-1;
            return 0;
        }
    }
    int ans=0;
    for(;i>=0;i--){
        Matrix ret=a*pw[i];
        if(ret.mat[0][0]-n<k){
            a=ret;
            if(ret.mat[0][0]-n<k) ans+=(1ll<<i);
        }
    }
    put(ans);
    return 0;
}

CF1067D Computer Game——矩阵优化与斜率优化，倍增二分思想

首先这个升级显然是吓唬你的，因为我可以一直选一个游戏 van，所以我只需要看 b_{i}p_{i} 最大就可以了。但是这里我们并不能考虑贪心，因为在时间短的情况下可能升级升不了，还是要 dp 的。

不难有 dp 方程如下，设 f(t) 表示还剩下 t 秒的最大期望，v 表示 b_{i}p_i 的最大值：

f(t+1)= \max_{i=1}^n \left\{ \underbrace{p_{i}(tv+a_i)}_{\text{升级成功}} + \underbrace{（1-p_{i})f_t}_{\text{升级失败}} \right\}

时间复杂度 O(nt)。

不难看出来可以斜率优化啊，但是我们要变下形式：

\begin{aligned} f(t+1) & = \max_{i=1}^n \left\{ p_{i}(tv+a_i) + (1-p_{i})f_t \right\} \\ & = p_{i}tv+p_{i}a_{i}+f_t-p_{i}f_{t} \\ & = p_{i}(tv-f_t)+p_{i}a_{i}+f_t \end{aligned}

因为 f_{t} 是已知的，所以这个就是一个显然的斜率优化式子，通过将 p_i 排序可以满足 k 单调，但是 x 呢？其实也是一样的：

\begin{aligned} tv-f_{t}& \ge (t-1)v-f_{t-1} \\ tv-f_{t}& \ge tv-v-f_{t-1} \\ f_{t-1}-f_{t} & \le v \end{aligned}

因为两个游戏之间获得的收益不可能比玩最大收益（最大的 b_{i}p_{i} 的游戏）还大，所以式子成立，x 单调不降。

故单调队列优化，时间复杂度 O(t+n).......t\le 10^{10}？

这个数据范围已经不行了，必须出矩阵优化....等会矩阵怎么斜率优化？

首先我们先把转移的矩阵搞出来，推啊推：

f_{i-1} & i-1 & 1 \end{bmatrix} \times \begin{bmatrix} (1-p_i) & 0 & 0\\ p_i v & 1 & 0\\ (p_i-a_i) & 1 & 1 \end{bmatrix} = \begin{bmatrix} f_{i} & i & 1 \end{bmatrix}

其实也不是很难推，有啥加啥，因为少个 1 直接加上去就行。

如果我们想找出来有哪些游戏是我们在斜率优化需要的，可以利用单调栈（不能用单调队列我们要存下来的）来记录我们斜率从那些点转移过来，现在问题在于如何确定什么时候从一个点转移到另一个点。

回忆一下这张图：

在斜率优化上，我们能用单调队列来做是因为对于每一个点上的斜率，它有一定转移的边界，在这之前是这个斜率，在之后就不是了。

说的好听矩阵怎么做？首先一个游戏的转移矩阵肯定不会变。问题在于我们怎么像单调队列优化一样找到所谓的边界呢？

首先单调队列不太行因为它不适用于矩阵这种玩意，那怎么办，矩阵这玩意也不能上斜率不单调三小将…………二分？

但其实维护凸壳的时候斜率函数单调递增，我们可以借助这个二分，找到顶点就可以了，其实二分也可以套在 k 与 x 同单调的地方，芝士没有那么优罢了

我们可以二分矩阵快速幂的幂，到哪个幂的时候转移是最优的！这样的时间复杂度是 O(n \log^2 t)，可以通过。

我们不妨快点，不难发现幂其实是一个倍增的形式，我们可以利用倍增的形式二分，首先预处理矩阵快速幂后的各个幂对应的矩阵，从大到小枚举倍增的幂，不断检查是否合法（即是否小于 t)，让后检查是否更优，直接赋值即可！时间复杂度 O(n \log t)，其中 \log t=33 可以通过。

代码如下，注意精度！！！！！

#include<bits/stdc++.h>
#define ll long long
#define double long double
using namespace std;
constexpr int MN=6e5+15;
constexpr double eps=1e-13;
struct Node{
    double k,b;
}ln[MN],cl[MN],s[MN];
ll n,t,top,tot,now;
double v;

struct Matrix{
    double mat[5][5];

    Matrix operator *(const Matrix &x)const{
        Matrix ret;
        memset(ret.mat,0,sizeof(ret.mat));
        for(int i=1;i<=3;i++){
            for(int j=1;j<=3;j++){
                for(int k=1;k<=3;k++){
                    ret.mat[i][j]+=mat[i][k]*x.mat[k][j];
                }
            }
        }
        return ret;
    }

}g[40],f;

bool cmp(Node x,Node y){
    if(fabs(x.k-y.k)<eps) return x.b<y.b;
    return x.k<y.k;
}

int ck(double x){
    if(fabs(x)<eps) return 0;
    return x>0?1:-1;
}

double gety(Node x,Node y){
    return (x.b-y.b)/(y.k-x.k);
}

int main(){
    cin>>n>>t;
    for(int i=1;i<=n;i++){
        double a,b,p;
        cin>>a>>b>>p;
        v=max(v,b*p);
        ln[i].k=p;
        ln[i].b=p*a;
    }
    sort(ln+1,ln+1+n,cmp);
    for(int i=1;i<=n;i++){
        // 先把那些相等斜率的全排了
        if(i == n || ck(ln[i].k - ln[i+1].k) != 0) cl[++tot]=ln[i];
    }
    for(int i=1;i<=tot;i++){
        // 单调栈处理转移的点
        while(top>1&&ck(gety(s[top],s[top-1])-gety(cl[i],s[top-1]))>=0) top--;
        s[++top]=cl[i];
    }
    // f 为 初始矩阵
    f.mat[1][3]=1;
    for(int i=1;i<=top;i++){
        double x=now*v-f.mat[1][1];
        while(i<top&&ck(x-gety(s[i],s[i+1]))>=0) i++;// 先把过时决策排了
        if(i<top) x=gety(s[i],s[i+1]);
        g[0].mat[1][2]=g[0].mat[1][3]=g[0].mat[2][3]=0;
        g[0].mat[2][2]=g[0].mat[3][2]=g[0].mat[3][3]=1;
        g[0].mat[1][1]=1-s[i].k,g[0].mat[2][1]=s[i].k*v,g[0].mat[3][1]=s[i].b;// 初始化矩阵
        for(int j=1;j<=35;j++) g[j]=g[j-1]*g[j-1];
        for(int j=35;j>=0;j--){
            ll np=now+(1ll<<j);
            if(np>=t) continue;
            // 如果决策更优或已经到头了即可转移
            if(i==top||ck(x-np*v+(f*g[j]).mat[1][1])>=0){
                f=f*g[j];
                now=np;
            }
        }
        f=f*g[0];
        now++;
        if(now==t) break;
    }
    cout<<fixed<<setprecision(10)<<f.mat[1][1];
    return 0;
}

4.后言

总结下来，矩阵优化的特征分为如下：

  1. 数据范围既有极大又有极小
  2. 具有线性递推特征，或者转移过来的状态数量偏小。
  3. DP 式子可以写成矩阵并且矩阵规模小
  4. 转移规模极大

用途，一个优化转移，一个优化边权较小图上的统计问题，基本步骤就是根据操作手册来就行。

转移过程中还是有一些 trick 的：

  1. 如果题目中有一些关键节点会改变转移矩阵，考虑到向量乘矩阵和矩阵乘矩阵不同复杂度，可以考虑强制中断，二进制分组处理。
  2. 如果转移矩阵仅由 0/1 构成，考虑 bitset 优化。
  3. 若图有边权，边权极小，可以考虑拆点
  4. 如果涉及到指数高精，想想性质或者 10 进制快速幂
  5. 矩阵快速幂也满足倍增的思想，可以考虑倍增二分

个人理解，矩阵优化的本质就是线性递推，其实回看所有拆点操作，都是将这些转移满足 i\rightarrow i+1 的形式，只有这样才能进行矩阵加速递推。

而对于操作手册是自己根据刷题经验总结出来的一些步骤，大家可以根据这几个步骤灵活调用。

完结撒花！