DP 优化方法大杂烩

Alex_Wei · 2021-07-29 17:23:26 · 个人记录

*标 的是推荐阅读的部分 / 做的题目**。

非常推荐在 cnblogs 中阅读。

1. 动态 DP（DDP）算法简介

动态动态规划。

以 P4719 为例讲一讲 ddp：

1.1. 树剖解法

如果没有修改操作，那么可以设计出 DP 方案 f_{i,0/1} 分别表示不选（0）/ 选（1）点 i 的最大权值，那么有 f_{i,0}=\sum_{x\in S_i}\max(f_{x,0},f_{x,1}),f_{i,1}=v_i+\sum_{x\in S_i}f_{i,0}。

如果加上修改操作，那么通常的一个想法是用矩阵表示转移方程再用数据结构维护。注意到一个节点可能有非常多的儿子，这样三次方级别的矩阵乘法就凉凉了。

此时我们就需要结合矩阵乘法 + 数据结构维护（通常是线段树）能够快速转移 DP 方程的优点，与树链剖分从根节点到任意一个节点所经过的轻重链个数为 \log 级别的优秀性质，设计一个新的算法。它应该满足从一个点的重儿子转移到这个点的方程可以用矩阵描述（想一想，为什么）。

首先对给出的树进行树链剖分，记 s_i 为 i 的重儿子。设 g_{i,0/1} 分别表示不选（0）/选（1）点 i 且只能从轻儿子中转移的最大权值，那么有 f_{i,0}=g_{i,0}+\max(f_{s_i,0},f_{s_i,1}),f_{i,1}=g_{i,1}+f_{s_i,0}。定义矩阵乘法 A\times B=C 为 C_{i,j}=\max(A_{i,k}+B_{k,j})。这样，我们就可以用 \begin{bmatrix}f_{s_i,0}&f_{s_i,1}\end{bmatrix}\times\begin{bmatrix}g_{i,0}&g_{i,1}\\g_{i,0}&-\infty\end{bmatrix} 得到 \begin{bmatrix}f_{i,0}&f_{i,1}\end{bmatrix}。f,g 需要在一开始预处理出来（其实 f 不用，但是预处理 g 依赖于求出 f）。记点 i 从其重儿子转移到该点的矩阵为 G_i，接下来考虑修改操作怎么办：

注意到一次修改操作只会改变 \log 级别的点的 G。这是因为只有在轻儿子被修改时需要修改 G，而在 i\to root（i 是被修改的点）的路径上只有 \log 条轻边。这样我们可以直接每次跳重链维护该重链顶部节点的父节点的 G。

那么怎么在更改一个点的轻儿子的 g 后维护该点的 G 呢（其实就是 g）？OI Wiki 上对于这个操作描述过于简略导致我苦思冥想了 15min 才搞懂。具体地，注意到 g 的转移方程相较于 f 其实只少了一个 x\neq s_i 的条件，而每个点对于该点父亲的 g 的贡献是以所有轻儿子的 f 的和的形式计算的，因此可以在修改前先计算原来的贡献，将其从 g_{top_i} 中减掉，再加上更新后的贡献即可。

注意到一个重链的尾端只有可能是叶子结点，而叶子结点没有 s_i，所以它的 f 值就等于 g。因此一个点的 f 值就等于它到它所在重链末端的矩阵乘积。记录一下 top=i 的重链的末端 ed_i，查询的时候直接求 \prod_{i=1}^{dfn_{ed_1}}G_i 即可。

注意点：因为 dfn 大的在底下，而又是从下往上转移，所以应该右区间乘左区间，查询的时候也要先跑右区间。调了 2147483647 h。

时间复杂度 \mathcal{O}(n\log^2 n)，代码见下方例题。常数有亿点点大。

1.2. LCT 解法

等我把 LCT 学会了再来填坑。

1.3. 例题

I. P4719 【模板】"动态 DP"&动态树分治

思路见最顶上的 DDP 算法简介，这里给出代码。

const int N=1e5+5;
const int inf=1e9+7;

struct matrix{
    int a,b,c,d;
    friend matrix operator * (matrix x,matrix y){
        matrix z;
        z.a=max(x.a+y.a,x.b+y.c);
        z.b=max(x.a+y.b,x.b+y.d);
        z.c=max(x.c+y.a,x.d+y.c);
        z.d=max(x.c+y.b,x.d+y.d);
        return z;
    }
}val[N<<2],G[N];
matrix prim,ans;

void build(int l,int r,int x){
    if(l==r)return val[x]=G[l],void();
    int m=l+r>>1;
    build(l,m,x<<1),build(m+1,r,x<<1|1);
    val[x]=val[x<<1|1]*val[x<<1]; 
} void modify(int l,int r,int p,int x){
    if(l==r)return val[x]=G[p],void();
    int m=l+r>>1;
    if(p<=m)modify(l,m,p,x<<1);
    else modify(m+1,r,p,x<<1|1);
    val[x]=val[x<<1|1]*val[x<<1];
} void query(int l,int r,int ql,int qr,int x){
    if(x==1)ans=prim;
    if(ql<=l&&r<=qr)return ans=ans*val[x],void();
    int m=l+r>>1;
    if(m<qr)query(m+1,r,ql,qr,x<<1|1);
    if(ql<=m)query(l,m,ql,qr,x<<1);
}

int n,q,dn,v[N],sz[N],son[N],fa[N];
int dfn[N],top[N],ed[N];
int f0[N],f1[N],g0[N],g1[N];
vector <int> e[N];
void gen(int x){
    int id=dfn[x];
    G[id].a=g0[x],G[id].b=g1[x];
    G[id].c=g0[x],G[id].d=-inf;
} void dfs1(int id,int f){
    fa[id]=f,sz[id]=1;
    for(int it:e[id])if(it!=f){
        dfs1(it,id);
        if(sz[son[id]]<sz[it])son[id]=it;
        sz[id]+=sz[it];
    }
} int dfs2(int id,int t){
    dfn[id]=++dn,top[id]=t,g1[id]=f1[id]=v[id];
    if(son[id]){
        ed[id]=dfs2(son[id],t);
        f1[id]+=f0[son[id]];
        f0[id]+=max(f0[son[id]],f1[son[id]]);
    } else ed[id]=id;
    for(int it:e[id])if(it!=fa[id]&&it!=son[id]){
        dfs2(it,it);
        g1[id]+=f0[it],g0[id]+=max(f0[it],f1[it]);
        f1[id]+=f0[it],f0[id]+=max(f0[it],f1[it]);
    } return ed[id];
} void modify(int x,int nv){
    g1[x]+=nv-v[x],v[x]=nv;
    while(top[x]!=1){
        int tp=top[x],ft=fa[tp];
        query(1,n,dfn[tp],dfn[ed[tp]],1);
        g1[ft]-=ans.a,g0[ft]-=max(ans.a,ans.b);
        gen(x),modify(1,n,dfn[x],1);
        query(1,n,dfn[tp],dfn[ed[tp]],1);
        g1[ft]+=ans.a,g0[ft]+=max(ans.a,ans.b),x=ft;
    } gen(x),modify(1,n,dfn[x],1);
} int query(){
    query(1,n,1,dfn[ed[1]],1);
    return max(ans.a,ans.b);
}

int main(){
    cin>>n>>q,prim.c=-inf,prim.d=-inf;
    for(int i=1;i<=n;i++)cin>>v[i];
    for(int i=1;i<n;i++){
        int u,v; cin>>u>>v;
        e[u].pb(v),e[v].pb(u);
    } dfs1(1,0),dfs2(1,1);
    for(int i=1;i<=n;i++)gen(i);
    build(1,n,1);
    for(int i=1;i<=q;i++){
        int x,y; cin>>x>>y;
        modify(x,y),cout<<query()<<endl;
    }
    return 0;
}

*II. CF750E New Year and Old Subsequence

题意简述：多次询问一个字符串的一个子串至少删去多少个字符后才能使得其不包含子序列 2016 而包含 2017。

动态 DP 板子题吧。感觉所有动态 DP 的题目都挺板的（flag）。

直接大力 DP 就好了：设 f_{i,0/1/2/3/4} 分别表示包含 \varnothing/2/20/201/2017 在 [1,i] 中最少需要删去多少个字符，那么有

需要注意的是如果 [] 中的条件不成立是不能转移的（而不是视为 0），应视为 \infty。因为转移只与前一项有关，故将上式改写成矩阵乘法形式即可。每次询问直接用 \begin{bmatrix}0&\infty&\infty&\infty&\infty\end{bmatrix} 去乘 [l,r] 的区间矩阵积（用线段树维护即可），得到的 f_{r,4} 即为所求（如果 f_{r,4}=\infty 则无解）。注意这里的矩阵乘法应定义为 C_{i,j}=\min(A_{i,k}+B_{k,j})。

事实上，\prod_{i=l}^r G_i 的第一行就是 f_r（G_i 是 f_{i-1}\to f_i 的转移矩阵），所以直接求出矩阵区间矩阵积即可，不过用上面的方法可以 f\times G 枚举第一行，矩阵乘法就是平方的；而 G\times G 是必须三方，**常数会更大。

/*
    Powered by C++11.
    Author : Alex_Wei.
*/

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

#define mem(x,v) memset(x,v,sizeof(x))

const int S=5;
const int N=2e5+5;
const int inf=0x3f3f3f3f;

struct matrix{
    int a[S][S];
    matrix friend operator * (matrix x,matrix y){
        matrix z; mem(z.a,0x3f);
        for(int i=0;i<S;i++)
            for(int j=0;j<S;j++)
                for(int k=0;k<S;k++)
                    z.a[i][j]=min(z.a[i][j],x.a[i][k]+y.a[k][j]);
        return z;
    }
    matrix friend operator / (matrix x,matrix y){
        matrix z; mem(z.a,0x3f);
        for(int i=0;i<S;i++)
            for(int j=0;j<S;j++)
                z.a[0][i]=min(z.a[0][i],x.a[0][j]+y.a[j][i]);
        return z;
    }
}val[N<<2],G[N],ans,pri;

void build(int l,int r,int x){
    if(l==r)return val[x]=G[l],void();
    int m=l+r>>1;
    build(l,m,x<<1),build(m+1,r,x<<1|1);
    val[x]=val[x<<1]*val[x<<1|1]; 
} void query(int l,int r,int ql,int qr,int x){
    if(ql<=l&&r<=qr)return ans=ans/val[x],void();
    int m=l+r>>1;
    if(ql<=m)query(l,m,ql,qr,x<<1);
    if(m<qr)query(m+1,r,ql,qr,x<<1|1);
}

int n,q;
char s;
int main(){
    cin>>n>>q;
    mem(pri.a,0x3f),pri.a[0][0]=0;
    for(int i=1;i<=n;i++){
        cin>>s,mem(G[i].a,0x3f);
        for(int j=0;j<S;j++)G[i].a[j][j]=0;
        if(s=='2')G[i].a[0][0]=1,G[i].a[0][1]=0;
        if(s=='0')G[i].a[1][1]=1,G[i].a[1][2]=0;
        if(s=='1')G[i].a[2][2]=1,G[i].a[2][3]=0;
        if(s=='7')G[i].a[3][3]=1,G[i].a[3][4]=0;
        if(s=='6')G[i].a[3][3]=1,G[i].a[4][4]=1;
    } build(1,n,1);
    for(int i=1;i<=q;i++){
        int l,r; cin>>l>>r;
        ans=pri,query(1,n,l,r,1);
        cout<<(ans.a[0][4]<=n?ans.a[0][4]:-1)<<endl;
    }
    return 0;
}

*III. P3781 [SDOI2017]切树游戏

题意简述：给出一个树 T，每个点有小于 2^7 的权值 v_i。多次查询 T 有多少个非空连通子树满足其所有点的权值异或和为 k。

先开个坑，到时候再来口胡。

想了一天感觉要用 FWT，然而我还不会。flag 倒了（悲）。看了眼题解，发现不在我的能力范围之内，先咕着再说。

2021.7.5 好的，今天学完 FWT 第一件事就是补这题。

首先设计 DP，设：

令 F_i(x),LF_i(x),H_i(x),LH_i(x) 分别为它们的生成函数（做过 FWT 后），那么有：

其中生成函数之间进行的是异或卷积。FWT，请。

综上，可以直接用矩阵描述转移：

树剖 + 线段树维护即可。修改时需要更新 LF 和 LH，LH 因为是和可以直接减，LF 可以维护乘以 0 的个数以及非零值之积。矩阵乘法的常数优化可见 基于变换合并的树上动态 DP 的链分治算法。视 n,q 同阶，则时间复杂度为 \mathcal{O}(nm\log^2 n)。

洛谷上树剖被卡了，LOJ 可过。

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

#define pb push_back
#define mcpy(x,y) memcpy(x,y,sizeof(y))
#define mem(x,v) memset(x,v,sizeof(x))

const int N=3e4+5;
const int M=1<<7;
const int mod=10007;

int n,m,q,v[N],inv[mod],tmp[M],V[N][M];
int F[N][M],LF[N][M],H[N][M],LH[N][M];
vector <int> e[N];

void print(string s,int *a){
    cout<<"output "<<s<<" : ";
    for(int i=0;i<m;i++)cout<<a[i]<<" ";
    cout<<endl;
}

// polynomial
void FWT(int *f,int x){
    for(int d=2,k=1;d<=m;d<<=1,k<<=1)
        for(int i=0;i<m;i+=d)
            for(int j=0;j<k;j++)
                f[i+j]=(f[i+j]+f[i+j+k])%mod,
                f[i+j+k]=(f[i+j]-f[i+j+k]-f[i+j+k]+mod+mod)%mod,
                f[i+j]=f[i+j]*x%mod,f[i+j+k]=f[i+j+k]*x%mod;
}
void init(int *a){for(int i=0;i<m;i++)a[i]=1;}
void mul(int *a,int *b){for(int i=0;i<m;i++)a[i]=a[i]*b[i]%mod;}
void mul(int *a,int *b,int *c){for(int i=0;i<m;i++)c[i]=a[i]*b[i]%mod;}
void mul_u(int *a,int *b){for(int i=0;i<m;i++)a[i]=a[i]*(b[i]+1)%mod;}
void add(int *a,int *b){for(int i=0;i<m;i++)a[i]=(a[i]+b[i])%mod;}
void sub(int *a,int *b){for(int i=0;i<m;i++)a[i]=(a[i]-b[i]+mod)%mod;}

int Z[N][M],X[N][M];
void upd_lf(int *g,int *x,int *z,int *b,bool tp){ // type = 1 : multiply ; type = 0 : divide
    for(int i=0;i<m;i++){
        int c=(b[i]+1)%mod;
        if(c){
            if(tp)x[i]=x[i]*c%mod;
            else x[i]=x[i]*inv[c]%mod;
        } else if(tp)z[i]++;
        else z[i]--;
        g[i]=z[i]?0:x[i];
    }
}

// tree partition
int dnum,sz[N],son[N],fa[N],ed[N],dep[N],dfn[N],top[N];
void dfs1(int id,int f,int d){
    fa[id]=f,dep[id]=d++,sz[id]=1;
    for(int it:e[id])
        if(it!=f){
            dfs1(it,id,d);
            if(sz[it]>sz[son[id]])son[id]=it;
            sz[id]+=sz[it];
        }
}
int dfs2(int id,int t){
    top[id]=t,dfn[id]=++dnum;
    mcpy(F[id],V[id]),init(LF[id]),init(X[id]);
    if(son[id]){
        ed[id]=dfs2(son[id],t);
        mul_u(F[id],F[son[id]]); // f
        add(H[id],H[son[id]]); // h
    } else ed[id]=id;

    for(int it:e[id])
        if(it!=fa[id]&&it!=son[id]){
            dfs2(it,it);
            upd_lf(LF[id],X[id],Z[id],F[it],1); // lf
            add(LH[id],H[it]); // lh
        }
    mul(F[id],LF[id]); // f
    add(H[id],LH[id]),add(H[id],F[id]); // h
    return ed[id];
}

// matrix
struct matrix{
    int a[M],b[M],c[M],d[M];
    friend matrix operator * (matrix x,matrix y){
        matrix z;
        mul(x.a,y.a,z.a);
        mul(x.a,y.b,z.b),add(z.b,x.b);
        mul(y.a,x.c,z.c),add(z.c,y.c);
        mul(x.c,y.b,z.d),add(z.d,x.d),add(z.d,y.d);
        return z;
    }
}G[N],val[N<<2],emp,ans;

void gen_mat(int x){
    int id=dfn[x];
    mul(LF[x],V[x],G[id].a);
    mcpy(G[id].b,G[id].a);
    mcpy(G[id].c,G[id].a);
    mcpy(G[id].d,G[id].a),add(G[id].d,LH[x]);
}

// segment tree
void push_up(int x){
    val[x]=val[x<<1|1]*val[x<<1];
}
void build(int l,int r,int x){
    if(l==r)return val[x]=G[l],void();
    int m=l+r>>1;
    build(l,m,x<<1),build(m+1,r,x<<1|1);
    push_up(x);
}
void modify(int l,int r,int p,int x){
    if(l==r)return val[x]=G[l],void();
    int m=l+r>>1;
    if(p<=m)modify(l,m,p,x<<1);
    else modify(m+1,r,p,x<<1|1);
    push_up(x);
}
void query(int l,int r,int ql,int qr,int x){
    if(ql<=l&&r<=qr){
        if(r==qr)ans=val[x];
        else ans=ans*val[x];
        return;
    }
    int m=l+r>>1;
    if(m<qr)query(m+1,r,ql,qr,x<<1|1);
    if(ql<=m)query(l,m,ql,qr,x<<1);
}

// modify
void modify(int x,int k){
    mem(V[x],0),V[x][v[x]=k]=1,FWT(V[x],1);
    while(top[x]!=1){
        int tp=top[x],ft=fa[tp];
        query(1,n,dfn[tp],dfn[ed[tp]],1);
        sub(LH[ft],ans.d),upd_lf(LF[ft],X[ft],Z[ft],ans.c,0);
        gen_mat(x),modify(1,n,dfn[x],1);
        query(1,n,dfn[tp],dfn[ed[tp]],1);
        add(LH[ft],ans.d),upd_lf(LF[ft],X[ft],Z[ft],ans.c,1),x=ft;
    } gen_mat(x),modify(1,n,dfn[x],1);
}

int main(){
    // init : inverse & unit
    inv[1]=1;
    for(int i=2;i<mod;i++)inv[i]=(-mod/i*inv[mod%i]%mod+mod)%mod;

    cin>>n>>m;
    int nm=1; while(nm<m)nm<<=1; m=nm;
    for(int i=1;i<=n;i++)cin>>v[i],V[i][v[i]]=1,FWT(V[i],1);
    for(int i=1,u,v;i<n;i++)cin>>u>>v,e[u].pb(v),e[v].pb(u);

    dfs1(1,1,1),dfs2(1,1);
    for(int i=1;i<=n;i++)gen_mat(i);
    build(1,n,1);

    cin>>q;
    for(int i=1,x,k;i<=q;i++){
        string s; cin>>s;
        if(s[0]=='C')cin>>x>>k,modify(x,k);
        else{
            query(1,n,1,dfn[ed[1]],1);
            FWT(ans.d,inv[2]);
            cin>>k,cout<<ans.d[k]<<endl;
        }
    }
    return 0;
}

2. 矩阵快速幂 / 矩阵乘法优化 DP

2.1. 算法简介

如果一个 DP，转移只与很少的状态有关，而转移的轮数很多，那很有可能就是矩阵快速幂 / 矩阵乘法优化 DP 了。

实际上动态 DP 就运用了这种思想。

2.2. 例题

*I. P3176 [HAOI2015]数字串拆分

题意简述：定义 f(s) 为将 s 拆分成若干个不大于 m 个数的方案数。给出数字字符串 t，求 g(t)=\sum f(x)，其中 x 为将 t 分割成若干个数后它们的和。例如 t=\texttt{12} 时，答案为 f(1+2)+f(12)。

在我的洛谷博客内查看。

注意到 m 的范围很小，只有 5，所以 f(s) 可以用矩阵快速幂 m^3\log s 计算（基本操作）。记转移矩阵为 G，那么答案即为 (\sum G^x)_{0,0}。我们记 d_i 为 g'(t[1:i])=\sum G^x（可以和上面 g(s) 的柿子对比一下，g 是数值，g' 是矩阵），则答案为 (d_n)_{0,0}。

不难求出转移方程 d_i=\sum_{j=0}^{i-1}d_j\times f(t[j+1:i])。因为 t[j+1:i] 会很大（达到 10^n），所以快速幂的时候需要对指数高精（矩阵没有费马小定理！）。因此，这样的时间复杂度为 \mathcal{O}(n^3m^3)，无法承受。

注意到 f(t)=G^t=\prod_i G^{10^{|t|-i}t_i}，那么我们可以预处理出形如 G^{c\times 10^k}\ (1\leq c<10,1\leq k\leq n) 的矩阵。然后记 D_{l,r} 为 f(t[l:r])，那么它可以由 D_{l+1,r}\times G^{10^{r-l}t_{l}} 转移过来。nm^3|\Sigma|+n^2m^3 预处理后可以 n^2m^3 DP（d_i=\sum_{j=0}^{i-1} d_jD_{j+1,i}）。可以通过本题。

/*
    Powered by C++11.
    Author : Alex_Wei.
*/

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

//#pragma GCC optimize(3)
//#define int long long

using ll = long long;
#define mem(x,v) memset(x,v,sizeof(x))

const int N=500+5;
const int M=5;
const int mod=998244353;

int n,m;
char s[N];
void add(int &x,ll y){
    x+=y%mod;
    if(x>mod)x-=mod;
}

struct mat{
    int a[M][M];
    friend mat operator * (mat x,mat y){
        mat z; mem(z.a,0);
        for(int i=0;i<m;i++)
            for(int j=0;j<m;j++)
                for(int k=0;k<m;k++)
                    add(z.a[i][j],1ll*x.a[i][k]*y.a[k][j]);
        return z;
    }
    friend mat operator + (mat x,mat y){
        mat z; mem(z.a,0);
        for(int i=0;i<m;i++)
            for(int j=0;j<m;j++)
                add(z.a[i][j],x.a[i][j]+y.a[i][j]);
        return z;
    }
}base,f[N],c[N][N],d[N][10];

mat ksm(mat a,int b){
    mat x; mem(x.a,0);
    for(int i=0;i<m;i++)x.a[i][i]=1;
    while(b){
        if(b&1)x=x*a;
        a=a*a,b>>=1;
    } return x;
}

int main(){
    scanf("%s%d",s+1,&m),n=strlen(s+1),f[0].a[0][0]=1;
    for(int i=0;i<m;i++)base.a[i][0]=1;
    for(int i=1;i<m;i++)base.a[i-1][i]=1;
    for(int i=0;i<10;i++){
        d[0][i]=ksm(base,i);
        for(int j=1;j<=n;j++)d[j][i]=ksm(d[j-1][i],10);
    } for(int r=n;r;r--)
        for(int l=r;l;l--)
            if(l==r)c[l][r]=d[0][s[r]-'0'];
            else c[l][r]=c[l+1][r]*d[r-l][s[l]-'0'];
    for(int i=1;i<=n;i++)
        for(int j=0;j<i;j++)
            f[i]=f[i]+f[j]*c[j+1][i];
    cout<<f[n].a[0][0]<<endl;
    return 0;
}

II. CF576D Flights for Regular Customers

将所有边按照 d_i 排序并依次加入。设初始状态为 A，那么有 A_{0,0}=1。设加入第 i 条边前邻接矩阵为 G_i，那么先 A\gets A\times G_i^{d_i-d_{i-1}}；然后加边，bfs 判断能否从可到达的点跑到 n。注意到邻接矩阵是 01 矩阵，于是可以用 bitset 优化。总时间复杂度 \frac{mn^3}{w}。

/*
    Powered by C++11.
    Author : Alex_Wei.
*/

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

#define mem(x,v) memset(x,v,sizeof(x))

const int N=150+5;
const int inf=0x3f3f3f3f;

int n,m,d[N],dis=inf;
struct mat{
    bitset <N> a[N];
    void clear(){
        for(int i=0;i<n;i++)
            for(int j=0;j<n;j++)
                a[i][j]=0;
    } friend mat operator * (mat x,mat y){
        mat z; z.clear();
        for(int i=0;i<n;i++)
            for(int j=0;j<n;j++)
                if(x.a[i][j])
                    z.a[i]|=y.a[j];
        return z;
    }
}base,ans;

struct edge{
    int u,v,w;
    bool operator < (const edge &x){
        return w<x.w;
    }
}e[N];

void ksm(mat &s,mat a,int b){
    while(b){
        if(b&1)s=s*a;
        a=a*a,b>>=1;
    }
}

bool vis[N];
int bfs(mat a){
    queue <pii> q; mem(vis,0);
    for(int i=0;i<n;i++)if(a.a[0][i])q.push({i,0}),vis[i]=1;
    while(!q.empty()){
        pii t=q.front(); q.pop();
        if(t.fi==n-1)return t.se;
        for(int i=0;i<n;i++)
            if(base.a[t.fi][i]&&!vis[i])
                q.push({i,t.se+1}),vis[i]=1;
    } return -1;
}

int main(){
    cin>>n>>m;
    for(int i=1;i<=m;i++)cin>>e[i].u>>e[i].v>>e[i].w;
    sort(e+1,e+m+1),ans.a[0][0]=1;
    if(e[1].w)cout<<"Impossible\n",exit(0);
    int step=0;
    for(int i=1;i<=m;i++){
        base.a[e[i].u-1][e[i].v-1]=1;
        int s=bfs(ans);
        if(s>=0)dis=min(dis,s+step);
        if(i<m)ksm(ans,base,e[i+1].w-e[i].w),step=e[i+1].w;
    } if(dis==inf)cout<<"Impossible\n";
    else cout<<dis<<endl;
    return 0;
}

III. P1707 刷题比赛

有点裸，不过可以用来熟悉熟悉矩阵快速幂。

/*
    Powered by C++11.
    Author : Alex_Wei.
*/

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

using ll = long long;

#define fi first
#define se second
#define mem(x,v) memset(x,v,sizeof(x))

const int N=32;
const int mod=1e9+7;

ll n,m,p,q,r,t,u,v,w,x,y,z;
void add(ll &x,ll y){x=(x+y)%m;}
ll mul(ll x,ll y){
    if(x<y)swap(x,y);
    ll ans=0,b=x;
    while(y){
        if(y&1)add(ans,b);
        y>>=1,add(b,b);
    } return ans;
}
struct mat{
    ll a[N][N];
    friend mat operator * (mat x,mat y){
        mat z; mem(z.a,0);
        for(int i=0;i<11;i++)
            for(int j=0;j<11;j++)
                for(int k=0;k<11;k++)
                    add(z.a[i][j],mul(x.a[i][k],y.a[k][j]));
        return z;
    }
}base,a;

int main(){
    cin>>n>>m>>p>>q>>r>>t>>u>>v>>w>>x>>y>>z;

    a.a[0][1]=a.a[0][3]=a.a[0][5]=a.a[0][6]=a.a[0][7]=a.a[0][10]=1;
    a.a[0][0]=a.a[0][2]=a.a[0][4]=3;
    a.a[0][8]=w,a.a[0][9]=z;

    base.a[0][0]=p,base.a[1][0]=q;
    base.a[2][0]=base.a[4][0]=1;
    base.a[6][0]=r,base.a[7][0]=t,base.a[10][0]=1;

    base.a[2][2]=u,base.a[3][2]=v;
    base.a[0][2]=base.a[4][2]=1;
    base.a[8][2]=1;

    base.a[4][4]=x,base.a[5][4]=y;
    base.a[0][4]=base.a[2][4]=1;
    base.a[9][4]=base.a[7][4]=1,base.a[10][4]=2;

    base.a[0][1]=base.a[2][3]=base.a[4][5]=base.a[10][10]=1;
    base.a[8][8]=w,base.a[9][9]=z;

    base.a[7][6]=2,base.a[6][6]=base.a[10][6]=1;
    base.a[7][7]=1,base.a[10][7]=1;

    n-=2;
    while(n){
        if(n&1)a=a*base;
        base=base*base,n>>=1;
    } printf("nodgd %lld\nCiocio %lld\nNicole %lld\n",a.a[0][0],a.a[0][2],a.a[0][4]);

    return 0;
}

*IV. P4569 [BJWC2011]禁忌

题解。

V. P5059 中国象棋

注意到每一列是独立的，且方案数为 fib_{n+3}-(n+2)：设 f_i 表示第 i 列（从 1 下标开始，共有 n+1 列）不放棋子且任意两个不放棋子的格子间隔不超过 2 的方案数，则有 f_0=f_1=1，f_i=f_{i-1}+f_{i-2}， 则 f_i=fib_{i+1}。总方案数为 f_{n+1}+f_n=f_{n+2}=fib_{n+3}，再减去全部不放的 1 种和只有一个放的 n+1 种即可。故答案为 (fib_{n+3}-n-2)^{n+1}。计算 fib 用矩阵快速幂即可。注意要用慢速乘。

时间复杂度 \mathcal{O}(\log^2 n)。

ll n, p;
void add(ll &x, ll y) {
    x = (x + y) % p;
}
ll mul(ll x, ll y) {
    ll s = 0;
    while(y) {
        if(y & 1) add(s, x);
        x = (x << 1) % p, y >>= 1;
    }
    return s;
}
ll ksm(ll a, ll b) {
    ll s = 1;
    while(b) {
        if(b & 1) s = mul(s, a);
        a = mul(a, a), b >>= 1;
    }
    return s;
}

struct Matrix {
    ll a[2][2];
    Matrix() {mem(a, 0, 2);}
    friend Matrix operator * (Matrix x, Matrix y) {
        Matrix z;
        for(int i = 0; i < 2; i++)
            for(int j = 0; j < 2; j++)
                for(int k = 0; k < 2; k++)
                    add(z.a[i][j], mul(x.a[i][k], y.a[k][j]));
        return z;
    }
} base, c, I;

Matrix ksm(Matrix a, ll b) {
    Matrix ans = I;
    while(b) {
        if(b & 1) ans = ans * a;
        a = a * a, b >>= 1;
    }
    return ans;
}

int main() {
    cin >> n >> p;
    I.a[0][0] = I.a[1][1] = 1;
    base.a[0][0] = base.a[1][0] = base.a[0][1] = 1;
    c.a[0][0] = 2, c.a[0][1] = 1;
    ll res = (c * ksm(base, n)).a[0][0] - n - 2;
    res = (res % p + p) % p;
    cout << ksm(res, n + 1) << endl;
    return 0;
}

VI. P1397 [NOI2013] 矩阵游戏

比较裸的矩阵加速。但是我们没有办法转成二进制（复杂度 \log^2 n 直接爆炸），那么直接十进制预处理即可。

时间复杂度 \mathcal{O}(\log n)。

const int N = 1e6 + 5;
const int mod = 1e9 + 7;

void read(short int *x, int &sz) {
    char s = getchar();
    while(!isdigit(s)) s = getchar();
    while(isdigit(s)) x[++sz] = s - '0', s = getchar();
    x[sz]--;
    for(int i = sz; ; i--)
        if(x[i] < 0) x[i] += 10, x[i - 1]--;
        else break;
    if(x[1] == 0) {
        for(int i = 1; i < sz; i++) x[i] = x[i + 1];
        sz--;
    }
}

struct Matrix {
    ll a, b, c, d;
    void clear() {a = b = c = d = 0;}
    void init() {a = d = 1, b = c = 0;}
    void print() {
        puts("print Matrix : ");
        cout << a << " " << b << endl;
        cout << c << " " << d << endl;
        puts("finished."); 
    }
    friend Matrix operator * (Matrix x, Matrix y) {
        Matrix z;
        z.a = (x.a * y.a + x.b * y.c) % mod;
        z.b = (x.a * y.b + x.b * y.d) % mod;
        z.c = (x.c * y.a + x.d * y.c) % mod;
        z.d = (x.c * y.b + x.d * y.d) % mod;
        return z;
    }
} col[15], res[15], row;

short int n[N], m[N];
int sn, sm;

void init(Matrix a, Matrix *pw) {
    pw[0].init(), pw[1] = a;
    for(int i = 2; i <= 10; i++) pw[i] = pw[i - 1] * a;
}
Matrix pow10(Matrix a) {
    Matrix b = a * a, c = b * b;
    return b * c * c;
}
Matrix ksm(Matrix *a, short int *x, int sz) {
    Matrix res; res.init();
    for(int i = 1; i <= sz; i++) {
        res = pow10(res);
        res = res * a[x[i]];
    } return res;
}

int main() {
    read(n, sn), read(m, sm);
    cin >> col[1].a >> col[1].c >> row.a >> row.c, col[1].d = row.d = 1;
    init(col[1], col);
    Matrix tmp = ksm(col, m, sm), line = tmp * row;
    init(line, res);
    Matrix ans = ksm(res, n, sn); ans = ans * tmp;
    cout << (ans.a + ans.c) % mod << endl;
    return 0;
}

3. 状压 DP

其实不能算是一种优化方法，不过也记在这了。

3.1. 例题

I. P1357 花园

注意到 m,k 非常的小而 n 非常的大，明示矩阵快速幂了。先考虑朴素 DP：状压一下，记 f_{i,S} 为位置 i 及其前 m-1 个位置所表示状态为 S 时的方案数，转移就枚举下一位是 C 还是 P，判断一下即可。注意是环形 DP，所以枚举一开始 m 个位置的状态 S，初始值 f_{m,S}=1，那么 ans\gets ans+f_{n+m,S}。

这个过程用矩阵快速幂优化一下即可，总时间复杂度 \mathcal{O}(2^m\log n(2^m)^3)=\mathcal{O}(16^m\log n)。

看了题解之后发现可以做到 8^m\log n，但是不是很会。

/*
    Powered by C++11.
    Author : Alex_Wei.
*/

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

using uint = unsigned int;
using ll = long long;

#define mem(x,v) memset(x,v,sizeof(x))
#define mcpy(x,y) memcpy(x,y,sizeof(y))

const int N=32;
const int mod=1e9+7;

ll n,m,k,sz,ans;
void add(ll &x,ll y){x=(x+y)%mod;}
int count(int x){int c=0; while(x)c+=x&1,x>>=1; return c;}
struct mat{
    ll a[N][N];
    friend mat operator * (mat x,mat y){
        mat z; mem(z.a,0);
        for(int i=0;i<sz;i++)
            for(int j=0;j<sz;j++)
                for(int k=0;k<sz;k++)
                    add(z.a[i][j],x.a[i][k]*y.a[k][j]);
        return z;
    }
}base,s;

mat ksm(mat s,mat a,ll b){
    while(b){
        if(b&1)s=s*a;
        a=a*a,b>>=1;
    } return s;
}

int main(){
    cin>>n>>m>>k,sz=(1<<m)-1;
    for(int i=0;i<1<<m;i++){
        if(count(i)>k)continue;
        ll a=(i<<1)&sz,b=((i<<1)+1)&sz;
        if(count(a)<=k)base.a[i][a]=1;
        if(count(b)<=k)base.a[i][b]=1;
    } for(int i=0;i<1<<m;i++){
        mem(s.a,0),s.a[0][i]=1;
        s=ksm(s,base,n),add(ans,s.a[0][i]);
    } cout<<ans<<endl;
    return 0;
}

*II. AT695 マス目

题解。

*III. ACM/ICPC Regional Aizu 2013 Hidden Tree

注意到一棵 balanced tree 上所有数都是最小数乘以 2 的幂，因此，对于不同的 p\times 2^r\ (2\nmid p) 我们单独处理：计算每个数中有多少个质因子 2，记做 r，状压 f_i 表示当前所有数的和为 i 时，最多的叶结点个数。按位置顺序枚举所有 2^r，并令 f_{k+2^r}\gets \max(f_{k+2^r},f_k+1)\ (2^r\mid k)，要有后面的条件是因为考虑加入 2^r 之后， 其后面比 2^r 小的部分就没法跨过 2^r 和前面比 2^r 小的部分互相抵消。记 c=\lfloor\log_2n\rfloor+\max r，其最大约为 8+9=17，则时间复杂度为 \mathcal{O}(Tn2^c)。

const int N=1e3+5;
const int M=1<<18;

int n,f[M];
vector <int> e[N];
void solve(){
    int ans=0;
    for(int i=1,r=1;i<=n;i++,r=1){
        int x=read();
        while(!(x&1))x>>=1,r<<=1;
        e[x].pb(r);
    }
    for(int i=1;i<N;i++)
        if(e[i].size()){
            int sum=0;
            for(int it:e[i]){
                for(int j=sum-sum%it;j>=0;j-=it)
                    f[j+it]=max(f[j+it],f[j]+1);
                sum+=it;
            }
            for(int i=1;i<=sum;i<<=1)ans=max(ans,f[i]);
            for(int i=1;i<=sum;i++)f[i]=-1e9-7;
            e[i].clear();
        }
    cout<<ans<<endl;
}
int main(){
    mem(f,0xcf,M),f[0]=0;
    while(n=read())solve();
    return 0;
}

4. 单调队列优化 DP

4.1. 单调队列优化多重背包

记物品个数为 n，背包容量为 V。每个物品体积为 V_i，价值为 val_i，最多放 c_i 个。f_{i,j} 为前 i 个物品体积之和为 j 的最大价值，那么暴力显然有

考虑每个 f_{i-1,j'} 对 f_{i,j}\ (j'\leq j,V_i|j-j') 的贡献：记 k=\frac{j-j'}{V_i}，那么 0\leq k\leq c_i，f_{i-1,j'}+k\times val_i\to f_{i,j}。注意到这个 k 与 j,j' 同时有关，那么在考虑相同的 i 不同的 j 的时候就比较麻烦。因为对于每个相差 V_i 的 j，相同的合法的 j' 的贡献不一样。

注意到当 j 增加 V_i 时，对于相同的 j'，其 k 值就会增加 1。一般遇到这样的情况，我们都会在 j' 处减去一定关于 j' 的贡献，并在 j 处增加一定关于 j 的贡献。考虑消除这样的影响，改写上式：记 j=qV_i+r\ (r<V_i)，则

注意到对于不同的 j，相同的 j'，q-k 的值是相同的（因为 j 增加 V_i 时，q,k 会同时增加 1），同时 q\times val_i 是一个只关于 j 的值，所以可以用单调队列优化前面的区间 \max。实际上，这就是将 j'\to j 的贡献减去了 \frac{j'}{V_i}，并将 j 自身的贡献增加了 \frac{j}{V_i}。可以看出这样并不改变 j'\to j 的贡献，因为 k=\frac{j-j'}{V_i}。

4.2. 例题

I. *P2254 [NOI2005] 瑰丽华尔兹

注意到求最长滑行距离就是求最少使用魔法的次数，那么设 f_{t,i,j} 表示 t 时间使钢琴在 (i,j) 的最少魔法使用次数。\mathcal{O}(nmt) 无法承受。

注意到一段时间内移动方向相同，那么设 f_{k,i,j} 表示在第 k 段结束后使钢琴在 (i,j) 的最少魔法使用次数。拿向右移动举例：设持续时间为 las=t_k-s_k+1，则

并且要满足 (i,p),(i,p+1),\cdots,(i,j) 之间没有障碍。

类似单调队列优化多重背包，将 len-(j-p) 拆成 p 和 len-j，这样就分别与 p,j 有关，可以使用单调队列优化。其余方向同理。时间复杂度为 \mathcal{O}(nmk)。

/*
    Powered by C++11.
    Author : Alex_Wei.
*/

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

#define mem(x,v) memset(x,v,sizeof(x))
#define mcpy(x,y) memcpy(x,y,sizeof(y))

const int N=200+5;
const int inf=1e9+7;
int n,m,x,y,k,s,t,len,dir,f[N][N];
int ans,hd,tl,d[N],g[N][N];
char c[N][N];

int main(){
    cin>>n>>m>>x>>y>>k,mem(f,0x3f),f[x][y]=0;
    for(int i=1;i<=n;i++)for(int j=1;j<=m;j++)cin>>c[i][j];
    while(k--){
        cin>>s>>t>>dir,len=t-s+1;
        if(dir==1)for(int j=1;j<=m;j++,hd=1,tl=0)for(int i=n;i;i--){
            if(c[i][j]=='x'){
                g[i][j]=inf,hd=1,tl=0;
                continue;
            } while(hd<=tl&&f[i][j]-i<=f[d[tl]][j]-d[tl])tl--; d[++tl]=i;
            while(hd<=tl&&d[hd]-len>i)hd++;
            g[i][j]=f[d[hd]][j]+len-d[hd]+i;
        } if(dir==2)for(int j=1;j<=m;j++,hd=1,tl=0)for(int i=1;i<=n;i++){
            if(c[i][j]=='x'){
                g[i][j]=inf,hd=1,tl=0;
                continue;
            } while(hd<=tl&&f[i][j]+i<=f[d[tl]][j]+d[tl])tl--; d[++tl]=i;
            while(hd<=tl&&d[hd]+len<i)hd++;
            g[i][j]=f[d[hd]][j]+len+d[hd]-i;
        } if(dir==3)for(int i=1;i<=n;i++,hd=1,tl=0)for(int j=m;j;j--){
            if(c[i][j]=='x'){
                g[i][j]=inf,hd=1,tl=0;
                continue;
            } while(hd<=tl&&f[i][j]-j<=f[i][d[tl]]-d[tl])tl--; d[++tl]=j;
            while(hd<=tl&&d[hd]-len>j)hd++;
            g[i][j]=f[i][d[hd]]+len-d[hd]+j;
        } if(dir==4)for(int i=1;i<=n;i++,hd=1,tl=0)for(int j=1;j<=m;j++){
            if(c[i][j]=='x'){
                g[i][j]=inf,hd=1,tl=0;
                continue;
            } while(hd<=tl&&f[i][j]+j<=f[i][d[tl]]+d[tl])tl--; d[++tl]=j;
            while(hd<=tl&&d[hd]+len<j)hd++;
            g[i][j]=f[i][d[hd]]+len+d[hd]-j;
        } mcpy(f,g);
    } for(int i=1;i<=n;i++)for(int j=1;j<=m;j++)ans=max(ans,t-f[i][j]);
    cout<<ans<<endl;
    return 0;
}

*5. wqs 二分

5.1. 算法详解

wqs 二分，全称忘情水王钦石二分，又称带权二分，斜率凸优化。它常见于限制选取物品个数的 DP 当中，有着很明显的标志，因此经常看起来比较套路。设 f(i) 表示最多选取 i 个物品时的答案，若 f 为凸函数则适用 wqs 二分。

引入：n 个数，最多分成 m 段，求每段的和的平方和的最小值。

• 算法 1：n^2m 暴力 DP。
• 算法 2：决策单调性优化 DP，可以做到 nm\log n。P.S. 关于决策单调性优化 DP，详见 command_block 的 blog，个人认为 ta 总结的很好。
• 算法 3：斜率优化，nm。

可以发现，上面三种算法的时间复杂度都和 m 有关，很讨厌，而一层 DP 的时间复杂度必定要和 n 有关，所以考虑对 m 做些手脚。注意到 f_{n,i}（下文简写为 f_i） 是一个关于 i 的凸函数，此时就可以使用 wqs 二分来降低复杂度了。

所以 wqs 二分究竟二分的是什么？是斜率。具体地，如果我们用一个斜率为 k 的直线去切这个凸包，会切到一个点，不妨设为 (p,f_p)。那么因为这个 f 是下凸的，所以 (p,f_p) 一定是在所有 (i,f_i) 中，斜率为 k 时在 y 轴截距最小的那一个，为 f_p-kp。如图所示。（图源，图中的 c 为相当于本文中的 -k)。

不难发现这就相当于每选取一段后将答案减去 -k，不限段数时取到的最小值 f_p 及其段数 p。这个可以在 n/n\log n 的时间复杂度内完成。

如果 p=m，那么说明我们已经找到了答案；如果 p<m，那么说明这个斜率还不够大（时刻注意 k 是负的）；如果 p>m，那么说明这个斜率还不够小（只是对于引入这一题，如果上凸则相反）。根据结果调整二分上下界即可。因为我们 DP 求得的是截距 d，二分求得的是斜率 k，所以最终答案为 d+km。

5.2. 细节

5.2.1. 数据类型

大部分题目因为 f_i 是整数，所以斜率 \dfrac{f_{i+1}-f_i}{i+1-i}=f_{i+1}-f_i 也是整数，因此 wqs 二分时上下边界和斜率 k 都是整数。

但是有些题目（在接下来的例题部分会有涉及）由于涉及小数运算，所以斜率并非一定是整数。故 wqs 二分时需要注意精度。

5.2.2. 三点共线

如果出现了三点共线的情况，这就可能导致我们无论如何也二分不到想要的 k。即斜率为 k 时 p<m，斜率为 k'=k+1（或 k-1，视题目类型而定）时 p'>m。

就拿引入中的例子来看：如图，如果求出了 k_2，那么无论是 B，C 还是 D，斜率为 k_2 时它们在 y 轴的截距 d 是相同的（三点共线），所以算出来的 d+k_2x_C 也一定相同，故不影响最终的答案。但是有一点需要注意（假设我们要找 k_2）：

• 如果在保证和最小的情况下，DP 出的是段数最大值，那么当 p<m 时应有 l\gets m+1，p>m 时应有 r\gets m。k=k_1 时截到的点是 B，因为 x_B<x_C，所以 k_1 一定不是我们要找的斜率，l\gets m+1；k=k_2 时截到的点是 D，因为 x_D>x_C，所以 k_2 可能是我们要找的斜率，r\gets m。所以一开始的 m 应为 \lfloor\dfrac{l+r}{2}\rfloor。

• 相反，如果 DP 出的是段数最小值，那么当 p<m 时应有 l\gets m，p>m 时应有 r\gets m-1。k=k_3 时截到的点是 D，因为 x_D>x_C，所以 k_3 一定不是我们要找的斜率，r\gets m-1；k=k_2 时截到的点是 B，因为 x_B<x_C，所以 k_2 可能是我们要找的斜率，l\gets m。所以一开始的 m 应为 \lfloor\dfrac{l+r}{2}\rfloor+1。

上述两种情况可以对比一下。因此，在保证答案最优的情况下，究竟是求物品个数最大值还是最小值，对二分的过程有很大的影响，这一点在写 wqs 二分时需要特别注意。

5.2.3. 答案的更新

如果边界值变为 m 则更新答案（直接给答案赋值而不是取 \min/\max）；相反，若边界值变为 m\pm1，则不更新答案，因为 m 没有被取到。

如果答案所对应的斜率没有被 DP 过怎么办？实际上，上下边界稍微扩大一点之后是不会出现这种情况的。读者自证不难。

5.2.4. Trick

可以用结构体将 DP 值与所选取的物品个数结合在一起，这样不仅方便更新 DP 值，还方便比较两个 DP 值大小，因为若 DP 值相同还需比较所选取的物品个数（重写运算符即可，具体可见下方例题）。

5.3. 参考文章

在此感谢这些博主。

• https://www.mina.moe/archives/6349/comment-page-1#comment-1923
• https://blog.csdn.net/a_forever_dream/article/details/105581221
• https://www.cnblogs.com/CreeperLKF/p/9045491.html

5.4. 例题

I. P2619 [国家集训队]Tree I

经典题。这实际上不是 wqs 二分优化 DP。

设 f_i 表示恰好取 i 条白边时的答案，不难发现 f_i 是下凸函数，因此可以使用 wqs 二分。

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

#define ll long long
#define mem(x,v) memset(x,v,sizeof(x))
#define mcpy(x,y) memcpy(x,y,sizeof(y))

const int N=5e4+5;

struct edge{
    int u,v,val,col;
    bool operator < (const edge &v) const{
        return val!=v.val?val<v.val:col<v.col;
    }
}e[N<<1];

int n,m,need,ans,f[N];
int find(int x){return f[x]==x?x:f[x]=find(f[x]);}

bool check(int mid){
    for(int i=1;i<=m;i++)if(e[i].col==0)e[i].val-=mid;
    int tot=0,val=0; sort(e+1,e+m+1);
    for(int i=0;i<n;i++)f[i]=i;
    for(int i=1;i<=m;i++){
        int u=find(e[i].u),v=find(e[i].v);
        if(u==v)continue;
        tot+=!e[i].col,val+=e[i].val,f[u]=v; 
    } 
    for(int i=1;i<=m;i++)if(e[i].col==0)e[i].val+=mid;
    if(tot>=need)return ans=val+need*mid,1;
    return 0;
}

int main(){
    cin>>n>>m>>need;
    for(int i=1;i<=m;i++)cin>>e[i].u>>e[i].v>>e[i].val>>e[i].col;
    int l=-102,r=102;
    while(l<r){
        int m=l+r>>1;
        if(check(m))r=m;
        else l=m+1;
    } cout<<ans<<endl;
    return 0;
}

II. CF739E Gosha is hunting - 洛谷

不难发现贡献函数满足四边形不等式，故具有决策单调性，直接决策单调性分治即可。时间复杂度 \mathcal{O}(nm\log n)。

ll n,m,a[N],f[N],g[N],s[N][N],p[N];
void solve(int l,int r,int ql,int qr){
    if(l>r)return;
    ll m=l+r>>1,ans=1e18,p=0,lim=min(m-1,(ll)qr);
    for(int i=ql;i<=lim;i++)
        if(ans>f[i]+s[i+1][m])
            p=i,ans=f[i]+s[i+1][m];
    g[m]=ans,solve(l,m-1,ql,p),solve(m+1,r,p,qr);
}
int main(){
    cin>>n>>m;
    while(n){
        for(int i=1;i<=n;i++)cin>>a[i],p[i]=p[i-1]+a[i];
        for(int i=1;i<=n;i++)for(int j=i+1;j<=n;j++)
            s[i][j]=s[i][j-1]+a[j]*(p[j-1]-p[i-1]);
        mem(f,0x3f,N),f[0]=0;
        for(int i=0;i<=m;i++)solve(1,n,0,n-1),cpy(f,g,N);
        cout<<f[n]<<endl,cin>>n>>m;
    }
    return 0;
}

7. 斜率优化

斜率优化是非常重要的一类 DP 优化方法。然而我鸽到了现在才学。

7.1. 算法介绍

如果一类最优化 DP 可以被写成 f_i=\min/\max f_j+cost_j+cost_i+F_iF_j，即一些只和 i,j 有关的项和一个和 i,j 都有关的项的和，那么一般就可以使用斜率优化。

具体地，把转移方程改写为 f_j+cost_j=F_iF_j+(f_i-cost_i)，设 y=f_j+cost_j，k=F_i，x=F_j，b=f_i-cost_i，那么它就是一个一次函数的表达式 y=kx+b。注意到 k 是固定的，而所有的 (x,y) 我们也已知，由于 f_i=b+cost_i，所以我们就是要找到这样的点 (x,y) 使得经过这个点的斜率为 k 的直线在 y 轴的截距最大 / 小（取决于是 \max 还是 \min）。

不妨设转移方程里面是 \min，这就相当于我们要动态维护一个下凸包，并每次在这个凸包上二分出斜率为 k 的直线切到的点，则这个点就是最优决策点。但一般情况下我们并不需要动态凸包，也不需要二分：

• 如果满足插入的点的纵坐标 x 有序，且每次查询的斜率 k 有序，那么可以单调队列维护凸包。具体方式可见 wqs 二分部分的例题 IX. Aliens。
• 如果只满足 x 有序，那么就需要单调栈维护凸包并在上面二分（例题 VII.）。
• 否则需要动态凸包 / 李超线段树。

那么斜率优化就讲完了。

注意特判斜率不存在的情况。

7.2. 优质文章

• https://www.cnblogs.com/terribleterrible/p/9669614.html。
• https://www.cnblogs.com/Xing-Ling/p/11210179.html。

7.3. 例题

wqs 二分部分的例题 IX. & X. & XI. & XII. 就需要使用斜率优化作为每次二分时内层的 DP。

I. P3195 [HNOI2008]玩具装箱

太经典了。设 s_i=\sum_{j=1}^i C_j+1 并将 L 增加 1，那么转移方程即为 f_i=\min_{j=0}^{i-1} f_j+(s_i-s_j-L)^2。

化简一下就是 \underline{f_j+s_j^2}_{\ y}=\underline{2(s_i-L)}_{\ k}\times\underline{s_j}_{\ x}+\underline{f_i-(s_i-L)^2}_{\ b}。

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

#define ll long long
#define ld long double

const int N=5e4+5;
const ld eps=1e-10;

ll n,L,hd,tl,d[N],f[N],s[N];
ll sq(ll x){return x*x;}
ll Y(int i){return f[i]+sq(s[i]);}
ll X(int i){return s[i];}
ld slope(int i,int j){return (ld)(Y(j)-Y(i))/(X(j)-X(i));}

int main(){
    cin>>n>>L,L++;
    for(int i=1;i<=n;i++)cin>>s[i],s[i]+=s[i-1]+1;
    d[hd=tl=1]=0;
    for(int i=1;i<=n;i++){
        while(hd<tl&&slope(d[hd],d[hd+1])+eps<=2*(s[i]-L))hd++;
        int j=d[hd]; f[i]=f[j]+sq(s[i]-s[j]-L);
        while(hd<tl&&slope(d[tl-1],d[tl])-eps>=slope(d[tl],i))tl--;
        d[++tl]=i;
    } cout<<f[n]<<endl;
    return 0;
}

II. P2120 [ZJOI2007]仓库建设

预处理 p_i 的前缀和 sp_i 和 1\sim i 的产品全部运送到位置 i 的代价 sv_i，则转移方程为 f_i=\min_{j=0}^{i-1}f_j+sv_i-sv_j+sp_j\times(x_i-x_j)。

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

#define ll long long
#define ld long double
#define gc getchar()

inline int read(){
    int x=0,sign=0; char s=gc;
    while(!isdigit(s))sign|=s=='-',s=gc;
    while(isdigit(s))x=(x<<1)+(x<<3)+(s-'0'),s=gc;
    return sign?-x:x;
}

const int N=1e6+5;
const ld eps=1e-10;

ll n,x[N],p[N],c[N],f[N],sp[N],sv[N];
ll X(int i){return sp[i];}
ll Y(int i){return f[i]-sv[i]+sp[i]*x[i];}
ll cal(int i,int j){return sv[j]-sv[i]-sp[i]*(x[j]-x[i])+c[j];}
ld slope(int i,int j){return (ld)(Y(j)-Y(i))/(X(j)-X(i));}

ll hd,tl,d[N];
int main(){
    cin>>n;
    for(int i=1;i<=n;i++){
        x[i]=read(),p[i]=read(),c[i]=read();
        sp[i]=sp[i-1]+p[i],sv[i]=sv[i-1]+sp[i-1]*(x[i]-x[i-1]);
    } hd=tl=1;
    for(int i=1;i<=n;i++){
        while(hd<tl&&slope(d[hd],d[hd+1])+eps<=x[i])hd++;
        f[i]=f[d[hd]]+cal(d[hd],i);
        while(hd<tl&&slope(d[tl-1],d[tl])-eps>=slope(d[tl],i))tl--;
        d[++tl]=i;
    } cout<<f[n]<<endl;
    return 0;
}

*III. P3648 [APIO2014]序列分割

一个非常关键的 \mathbf{Observation} 是切的顺序无关，证明的话很简单，假设最终切出来 k+1 个块的元素和为 c_1,c_2,\cdots,c_{k+1}，那么答案即为 \sum_{i=1}^{k+1}\sum_{j=i+1}^{k+1}c_ic_j。因此可以愉快 DP 了：设 s_i 为 a_i 的前缀和，f_{i,p} 为前 i 个数分 p 个块的最大值，那么有 f_{i,p}=\max_{j=0}^{i-1}f_{j,p-1}+s_j\times (s_i-s_j)。注意记录每个 f_{i,p} 的最优决策点，方便输出方案，并特判 s_i=s_j 的情况（因为 a_i 可能等于 0）。时间复杂度 \mathcal{O}(nk)，使用 long double 会被卡常。

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

#define ll long long
#define ld double

const int N=1e5+5;
const int K=200+5;
const ld eps=1e-9;
const ld inf=1e9;

ll n,k,hd,tl,d[N],a[N],s[N],f[N],g[N];
int tr[N][K];
ld sl(int i,int j,int k){
    ld dx=s[j]-s[i],dy=f[j]-s[j]*s[j]-f[i]+s[i]*s[i];
    return fabs(dx)<eps?-inf:dy/dx;
} void print(int p,int t){
    if(!t)return;
    print(tr[p][t-1],t-1),cout<<p<<" ";
}
int main(){
    cin>>n>>k;
    for(int i=1;i<=n;i++)cin>>a[i],s[i]=s[i-1]+a[i];
    for(int i=1;i<=k;i++){
        d[hd=tl=1]=i-1;
        for(int j=i;j<=n;j++){
            while(hd<tl&&sl(d[hd],d[hd+1],i)+eps>=-s[j])hd++;
            int fr=d[hd];
            tr[j][i]=fr,g[j]=f[fr]+s[fr]*(s[j]-s[fr]);
            while(hd<tl&&sl(d[tl-1],d[tl],i)+eps<=sl(d[tl],j,i))tl--;
            d[++tl]=j;
        }
        memcpy(f,g,sizeof(g));
    } cout<<f[n]<<endl,print(tr[n][k],k);
    return 0;
}

IV. P6047 丝之割

具体地，如果两条弦 (u_i,v_i) 和 (u_j,v_j) 满足 u_i\leq u_j 且 v_i\geq v_j，那么 j 显然是无用的，因为割掉 i 的同时也一定能割掉 j。因此有用的弦一定满足 u_i,v_i 同时单调递增。像极了 “IOI2016 aliens”。

设 p_i 为 a_i 的前缀最小值，q_i 为 b_i 的后缀最小值，f_i 为割掉前 i 条弦的最小代价，那么有 f_i=\min_{j=0}^{i-1}f_j+a_{u_{j+1}-1}\times b_{v_i-1}。接下来斜率优化即可。时间复杂度 n+m\log m，瓶颈在排序，可以使用桶排优化到线性。

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

#define ll long long
#define ld double

const int N=3e5+5;
const ld eps=1e-9;

ll n,m,cnt,p[N],q[N],f[N];
struct str{
    int u,v;
    bool operator < (const str &x) const{
        return u==x.u?v>x.v:u<x.u;
    }
}c[N];

int hd,tl,d[N];
ld sl(int i,int j){
    ld dx=-p[c[j+1].u-1]+p[c[i+1].u-1];
    if(fabs(dx)<eps)return f[j]-f[i]<0?-1e9:1e9;
    else return (ld)(f[j]-f[i])/dx;
}

int main(){
    cin>>n>>m;
    for(int i=1;i<=n;i++)cin>>p[i];
    for(int i=1;i<=n;i++)cin>>q[i];
    for(int i=1;i<=m;i++)cin>>c[i].u>>c[i].v;
    for(int i=2;i<=n;i++)p[i]=min(p[i],p[i-1]);
    for(int i=n-1;i;i--)q[i]=min(q[i],q[i+1]);
    sort(c+1,c+m+1);
    for(int i=1,r=0;i<=m;i++)if(c[i].v>r)c[++cnt]=c[i],r=c[i].v;
    m=cnt,hd=tl=1;
    for(int i=1;i<=m;i++){
        while(hd<tl&&sl(d[hd],d[hd+1])+eps<=q[c[i].v+1])hd++;
        f[i]=f[d[hd]]+p[c[d[hd]+1].u-1]*q[c[i].v+1];
        while(hd<tl&&sl(d[tl-1],d[tl])-eps>=sl(d[tl],i))tl--;
        d[++tl]=i;
    } cout<<f[m]<<endl;
    return 0;
}

V. CF311B Cats Transport

经典题。首先将所有 t_i 减去 \sum_{j=2}^{H_i}D_j 求出如果要恰好接到第 i 只猫要从哪一时刻出发，然后设 s_i 为 t_i 的前缀和，f_{i,j} 表示用 j 个人接到前 i 只猫的最小等待时间之和，那么有 f_{i,j}=\min_{k=0}^{i-1}f_{k,j-1}+(i-k)t_i-(s_i-s_k)，斜率优化即可。

注意出发时间可以为负数，即不需要特判 t_i<0 的情况。

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

#define ll long long
#define ld long double
#define pb push_back
#define pii pair <int,int>
#define fi first
#define se second
#define gc getchar()
#define mcpy(x,y) memcpy(x,y,sizeof(y))
#define mem(x,v) memset(x,v,sizeof(x))

inline int read(){
    int x=0; char s=gc;
    while(!isdigit(s))s=gc;
    while(isdigit(s))x=x*10+s-'0',s=gc;
    return x;
}

const int N=1e5+5;
const ld eps=1e-9;

int n,m,p,hd,tl,d[N];
ll t[N],s[N],f[N],g[N];

ld sl(int i,int j){return (ld)(g[j]+s[j]-g[i]-s[i])/(j-i);}
int main(){
    cin>>n>>m>>p;
    for(int i=2;i<=n;i++)cin>>d[i],d[i]+=d[i-1];
    for(int i=1,h;i<=m;i++)cin>>h>>t[i],t[i]-=d[h];
    sort(t+1,t+m+1);
    for(int i=1;i<=m;i++)s[i]=s[i-1]+t[i];
    mem(g,0x3f),g[0]=0;
    for(int z=1;z<=p;z++){
        d[hd=tl=1]=0;
        for(int i=1;i<=m;i++){
            while(hd<tl&&sl(d[hd],d[hd+1])+eps<=t[i])hd++;
            int j=d[hd]; f[i]=g[j]+(i-j)*t[i]-s[i]+s[j];
            while(hd<tl&&sl(d[tl-1],d[tl])-eps>=sl(d[tl],i))tl--;
            d[++tl]=i;
        } mcpy(g,f);
    } cout<<f[m]<<endl;
    return 0;
}

VI. P2365 任务安排

经典题。发现一组任务的完成时间依赖于分成的组数，很难受，那么使用 代价提前计算 的技巧，即将启动时间 s 对 j+1\sim n 的所有任务的代价计算在内，而不仅仅是对于 j+1\sim i 的任务。

因此有转移方程 f_i=\min_{j=0}^{i-1}f_j+s(sf_n-sf_j)+st_i(sf_i-sf_j)。带个 s 是前缀和。

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

#define ll long long
#define ld long double
#define pb push_back
#define pii pair <int,int>
#define fi first
#define se second
#define gc getchar()
#define mcpy(x,y) memcpy(x,y,sizeof(y))
#define mem(x,v) memset(x,v,sizeof(x))

inline int read(){
    int x=0; char s=gc;
    while(!isdigit(s))s=gc;
    while(isdigit(s))x=x*10+s-'0',s=gc;
    return x;
}

const int N=5e3+5;
const ld eps=1e-10;

ll n,s,f[N],t[N],g[N],hd,tl,d[N];
ld sl(int i,int j){return (ld)(g[j]-s*f[j]-g[i]+s*f[i])/(f[j]-f[i]);}
int main(){
    cin>>n>>s;
    for(int i=1;i<=n;i++)cin>>t[i]>>f[i],f[i]+=f[i-1],t[i]+=t[i-1];
    hd=tl=1;
    for(int i=1;i<=n;i++){
        while(hd<tl&&sl(d[hd],d[hd+1])+eps<=t[i])hd++;
        int j=d[hd]; g[i]=g[j]+s*(f[n]-f[j])+t[i]*(f[i]-f[j]);
        while(hd<tl&&sl(d[tl-1],d[tl])-eps>=sl(d[tl],i))tl--;
        d[++tl]=i;
    } cout<<g[n]<<endl;
    return 0;
}

*VII. P5785 [SDOI2012]任务安排

和上题差不多。注意到 T_i 可能是负数，所以查询的斜率不一定单调增，但是插入的点的纵坐标一定单调不降（注意特判纵坐标相同的情况，因为 C_i 可能等于 0），所以要单调栈维护凸包，查询时二分斜率。

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

#define ll long long
#define ld long double
#define pb push_back
#define pii pair <int,int>
#define fi first
#define se second
#define gc getchar()
#define mcpy(x,y) memcpy(x,y,sizeof(y))
#define mem(x,v) memset(x,v,sizeof(x))

inline int read(){
    int x=0; char s=gc;
    while(!isdigit(s))s=gc;
    while(isdigit(s))x=x*10+s-'0',s=gc;
    return x;
}

const int N=3e5+5;
const ld eps=1e-10;
const ld inf=1e18;

ll n,s,f[N],t[N],g[N],hd,tl,d[N];
ld sl(int i,int j){
    ll dy=g[j]-s*f[j]-g[i]+s*f[i];
    if(f[i]==f[j])return dy<0?-inf:inf;
    return (ld)dy/(f[j]-f[i]);
}
int main(){
    cin>>n>>s;
    for(int i=1;i<=n;i++)cin>>t[i]>>f[i],f[i]+=f[i-1],t[i]+=t[i-1];
    hd=tl=1;
    for(int i=1;i<=n;i++){
        int l=hd,r=tl;
        while(l<r){
            int mid=l+r>>1;
            if(sl(d[mid],d[mid+1])+eps<=t[i])l=mid+1;
            else r=mid;
        } int j=d[l];
        g[i]=g[j]+s*(f[n]-f[j])+t[i]*(f[i]-f[j]);
        while(hd<tl&&sl(d[tl-1],d[tl])-eps>=sl(d[tl],i))tl--;
        d[++tl]=i;
    } cout<<g[n]<<endl;
    return 0;
}

*VIII. P2305 [NOI2014] 购票

首先将 s_i 前缀和处理，那么 g_i=\min_{j\in \mathrm{anc}(i)\land s_i-s_j\leq l_i}f_j+(s_i-s_j)p_i+q_i。斜率优化显然。

注意到每次查询一段后缀的凸包，而且是树上询问需要撤销，那么用线段树 / BIT 套可撤销单调栈维护凸包即可。由于 BIT 查询的是一段前缀，所以需要翻转一下。可以预处理 BIT 上每个节点维护的单调栈所需要的数组大小，这样就不需使用 vector 动态开空间（会 MLE）。时间复杂度 \mathcal{O}(n\log^2 n)。

类似题目：CEOI2009 Harbingers。本题见博客 “线段树的高级应用” 李超树部分例题 IV.

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

#define ll long long
#define ld long double

const int N=2e5+5;
ll n,t,fa[N],s[N],c[N],f[N];
ll p[N],q[N],l[N];
vector <int> e[N];
int get(ll x){return lower_bound(c+1,c+n+1,x)-c;}
ld sl(int i,int j){return (ld)(f[j]-f[i])/(s[j]-s[i]);}

int sz[N],ts[N],td[N],st[N<<4];
pair <int,int> d[N<<4];
void push(int id,int p){
    int b=sz[id]+ts[id];
    while(ts[id]>1&&sl(st[b-2],st[b-1])>sl(st[b-1],p))
        d[sz[id]+td[id]++]={p,st[b-1]},ts[id]--,b--;
    st[b++]=p,ts[id]++;
}
void del(int id,int p){
    ts[id]--; int b=sz[id]+td[id];
    while(td[id]&&d[b-1].first==p)st[sz[id]+ts[id]++]=d[b-1].second,td[id]--,b--;
}
int query(int id,ld x){
    if(!ts[id])return -1;
    int l=sz[id],r=sz[id]+ts[id]-1;
    while(l<r){
        int m=l+r>>1;
        if(sl(st[m],st[m+1])<=x)l=m+1;
        else r=m;
    } return st[l];
}

void dfs(int id){
    int x=n-get(s[id]-l[id])+1,y=n-get(s[id])+1,z=y;
    if(id>1)f[id]=1e18;
    while(x){
        int pos=query(x,p[id]);
        if(pos!=-1)f[id]=min(f[id],f[pos]+(s[id]-s[pos])*p[id]+q[id]);
        x-=x&-x;
    } while(y<=n)push(y,id),y+=y&-y;
    for(int it:e[id])dfs(it);
    while(z<=n)del(z,id),z+=z&-z;
}
int main(){
    cin>>n>>t;
    for(int i=1;i<=n;i++)for(int j=i;j<=n;j+=j&-j)sz[j+1]++;
    for(int i=3;i<=n;i++)sz[i]+=sz[i-1];
    for(int i=2;i<=n;i++){
        scanf("%lld%lld%lld%lld%lld",&fa[i],&s[i],&p[i],&q[i],&l[i]);
        c[i]=s[i]+=s[fa[i]],e[fa[i]].push_back(i);
    } sort(c+1,c+n+1),dfs(1);
    for(int i=2;i<=n;i++)printf("%lld\n",f[i]);
    return 0;
}

IX. P2900 [USACO08MAR]Land Acquisition G

显然如果两块土地 i,j 满足 w_i\leq w_j 且 l_i\leq l_j 那么 i 完全无用。因此将所有土地按照 w_i 从大到小排序，筛选出所有有用的土地，那么答案即为 f_i=\min_{0\leq j<i}f_j+w_{j+1}\times l_i。斜率优化即可。

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

#define ll long long
#define ld long double

const int N=5e4+5;

ll n,q,hd,tl,d[N],f[N];
struct seq{
    ll w,l;
    bool operator < (const seq &v) const{
        return w>v.w;
    }
}a[N],b[N];
ld sl(int i,int j){return (ld)(f[j]-f[i])/(b[i+1].w-b[j+1].w);}

int main(){
    cin>>n;
    for(int i=1;i<=n;i++)cin>>a[i].w>>a[i].l;
    sort(a+1,a+n+1);
    for(int i=1,r=0;i<=n;i++)if(a[i].l>r)r=a[i].l,b[++q]=a[i];
    for(int i=1;i<=q;i++){
        while(hd<tl&&sl(d[hd],d[hd+1])<=b[i].l)hd++;
        f[i]=f[d[hd]]+b[d[hd]+1].w*b[i].l;
        while(hd<tl&&sl(d[tl-1],d[tl])>=sl(d[tl],i))tl--;
        d[++tl]=i;
    } cout<<f[q]<<endl;
    return 0;
}

*X. P3571 [POI2014]SUP-Supercomputer

题意简述：一棵以 1 为根的树。q 次询问，每次给出 k，求至少要多少次同时访问不超过 k 个父节点已经被访问过的节点，才能访问完整棵树。根节点无限制。

sweet tea.

对于每一个 k，一定存在 i 使得深度不大于 i 的节点用 i 次访问，且深度大于 i 的节点每次都能访问 k 个（除了最后一次）。记 s_i 表示深度不小于 i 的节点个数，答案即为 \max_{i=1}^d\left(i+\lceil\dfrac{s_{i+1}}k\rceil\right)，其中 d 是最大深度。

若 i 是最优决策，那么对于任意一个 j\neq i，有 i+\lceil\dfrac{s_{i+1}}k\rceil\geq j+ \lceil\dfrac{s_{j+1}}k\rceil。略作变形得到 i-j \geq \lceil\dfrac{s_{j+1}-s_{i+1}}k\rceil。令横坐标为深度，纵坐标为 s_{x+1}，再写成斜率的形式，即当 j<i 时，\dfrac{s_{i+1}-s_{j+1}}{i-j}\geq -k，当 j>i 时，\dfrac{s_{i+1}-s_{j+1}}{i-j}\leq -k。不难看出这是一个上凸包的形式，即斜率递减。

具体地，我们对 (i,s_i) 建出上凸包，然后当 k 递增时，-k 递减，顶点会向横坐标大的方向移动，用指针维护即可。时间复杂度 \mathcal{O}(n)。

经过卡常拿到了最优解。

const int N=1e6+5;

int n,q,mxd,mxq,dep[N],f[N],qu[N];
int d[N],hd=1,tl;
ll s[N];

int main(){
    cin>>n>>q,dep[1]=s[1]=1;
    for(int i=1;i<=q;i++)mxq=max(mxq,qu[i]=read());
    for(int i=2,a;i<=n;i++)mxd=max(mxd,dep[i]=dep[read()]+1),s[dep[i]]++;
    for(int i=mxd-1;i;i--)s[i]+=s[i+1];
    for(int i=0;i<=mxd;i++){
        while(hd<tl&&(s[d[tl]+1]-s[d[tl-1]+1])*(i-d[tl])<=(s[i+1]-s[d[tl]+1])*(d[tl]-d[tl-1]))tl--;
        d[++tl]=i;
    }
    for(int i=1;i<=mxq;i++){
        while(hd<tl&&s[d[hd+1]+1]-s[d[hd]+1]>-i*(d[hd+1]-d[hd]))hd++;
        f[i]=d[hd]+(d[hd]==mxd?0:((s[d[hd]+1]-1)/i+1));
    }
    for(int i=1;i<=q;i++)print(f[qu[i]]),pc(' ');
    return flush(),0;
}

XI. P2497 [SDOI2012]基站建设

首先求出 j 对于 i 的贡献：根据勾股定理有 (x_i-x_j)^2+(r_{2,i}-r_{1,j})^2=(r_{2,i}+r_{1,j})^2。略作变形得到 4r_{2,i}r_{1,j}=(x_i-x_j)^2，故 \sqrt{r_{2,i}}=\dfrac{x_i-x_j}{2\sqrt{r_{1,j}}}。

显然的斜率优化，f_j-\dfrac {x_j}{2\sqrt{r_{1,j}}}\ (y_j)=x_i\times \left(-\dfrac{1}{2\sqrt{r_{1,j}}}\right)+(f_i-v_i)。注意到插入的 -\dfrac 1 {2\sqrt {r_{1,j}}} 并不单调，所以将其相反数看做斜率 k，y_j 看做截距插入动态开点李超线段树即可。时间复杂度 \mathcal{O}(n\log x)。

const int N = 5e5 + 5;
const ld inf = 1e18;

ll  m, x[N];
int n, R, node, mx[N << 5], ls[N << 5], rs[N << 5];
ld r[N], v[N], f[N], k[N], b[N];
ld get(int id, ld p) {return k[id] * p + b[id];}
void modify(ll l, ll r, int &x, int v) {
    if(!x) x = ++node;
    if(l == r) {
        if(get(mx[x], l) > get(v, l)) mx[x] = v;
        return;
    }
    ll m = l + r >> 1;
    if(get(mx[x], m) > get(v, m)) swap(v, mx[x]);
    if(get(mx[x], l) > get(v, l)) modify(l, m, ls[x], v);
    else modify(m + 1, r, rs[x], v);
}
ld query(ll l, ll r, int x, ll p){
    int m = l + r >> 1; ld ans = get(mx[x], p);
    if(l == r || !x) return ans;
    return min(ans, p <= m ? query(l, m, ls[x], p) : query(m + 1, r, rs[x], p));
}

int main() {
    cin >> n >> m, b[0] = inf;
    for(int i = 1; i <= n; i++) cin >> x[i] >> r[i] >> v[i];
    for(int i = 1; i <= n; i++) {
        if(r[i] == 0) {
            f[i] = inf;
            continue;
        }
        if(i == 1) f[i] = v[i];
        else f[i] = query(x[1], x[n], R, x[i]) + v[i];
        k[i] = 0.5 / sqrt(r[i]), b[i] = f[i] - x[i] / (2 * sqrt(r[i]));
        modify(x[1], x[n], R, i);
    }
    ld ans = inf;
    for(int i = 1; i <= n; i++)
        if(x[i] + r[i] >= m)
            ans = min(ans, f[i]);
    printf("%.3LF\n", ans);
    return 0;
}

XII.P4655 [CEOI2017]Building Bridges

有一个显然的 n^2 DP：f_i=\min_{\\j=1}^{i-1}f_j+\left(h_i-h_j\right)^2+\sum_{k=j+1}^{i-1}w_k，后面那一坨前缀和优化一下，再把平方拆开，移个项，得

显然的斜率优化形式，注意斜率和下标都不单调，故使用李超线段树维护即可。

const int N = 1e5 + 5;
const int H = 1e6;
const ll inf = 1e18;

ll n, f[N], h[N], w[N];

// Li Chao Segment Tree
int R, node, ls[N << 5], rs[N << 5], mi[N << 5];
ll k[N], b[N];
ll get(int x, int id) {return k[id] * x + b[id];}
void modify(int l, int r, int &x, int v) {
    if(!x) x = ++node;
    if(l == r) {
        if(get(l, v) < get(l, mi[x])) mi[x] = v;
        return;
    } int m = l + r >> 1;
    if(get(m, v) < get(m, mi[x])) swap(mi[x], v);
    if(get(l, v) < get(l, mi[x])) modify(l, m, ls[x], v);
    else if(get(r, v) < get(r, mi[x])) modify(m + 1, r, rs[x], v);
}
ll query(int l, int r, int p, int x) {
    if(l == r || !x) return get(p, mi[x]);
    ll m = l + r >> 1, ans = get(p, mi[x]);
    if(p <= m) return min(query(l, m, p, ls[x]), ans);
    return min(query(m + 1, r, p, rs[x]), ans);
}

int main() {
    cin >> n, mem(f, 63, N), f[1] = 0, b[0] = inf;
    for(int i = 1; i <= n; i++) cin >> h[i];
    for(int i = 1; i <= n; i++) cin >> w[i], w[i] += w[i - 1];
    for(int i = 2; i <= n; i++) {
        int j = i - 1;
        k[j] = -2 * h[j], b[j] = f[j] + h[j] * h[j] - w[j];
        modify(-H, H, R, j);
        f[i] = query(-H, H, h[i], R) + h[i] * h[i] + w[j];
    }
    cout << f[n] << endl;
    return 0;
}

XIII. P6302 [NOI2019] 回家路线 加强版

经典斜率优化。注意到如果按照站点 DP 则还需记录时间这一维，无法接受，不如按照列车 DP：设 f_i 表示搭乘第 i 号列车的最小代价。把柿子写出来：

时间造成的烦躁值可以在统计答案时算上，也可以算在 DP 方程里面。

注意到时间所产生的偏序关系为 DP 定下了一个顺序：按照时间 DP，将每号列车拆成出发和到达两个事件并按时间排序。注意因为可以刚下车就上车，即同一时刻从到达转移到出发，所以对于时刻相同的两个事件，到达应该排在出发之前。

由于时刻递增即斜率和插入点的横坐标有序，我们只需要对每个站点维护一个用单调队列实时维护的下凸壳即可。时间复杂度 \mathcal{O}(m\log m)，复杂度瓶颈在于排序。若使用桶排可做到 \mathcal{O}(m)。

时刻注意对于横坐标相同，斜率不存在的情况的处理。

const int N = 1e5 + 5;
const int M = 1e6 + 5;
const ld eps = 1e-9;
const ll inf = 1e15;

// type = 1 : depart
// type = 2 : arrive

struct Event {
    ll id, type, sta, time;
    bool operator < (const Event &v) const {
        return time != v.time ? time < v.time : type > v.type;
    }
} tr[M << 1];
deque <pll> conv[N];
ll n, m, A, B, C, ans = inf;
ll x[M], y[M], p[M], q[M], res[M];

ld Slope(pll a, pll b) {
    if(a.fi == b.fi) return b.se > a.se ? inf : -inf;
    return (ld)(b.se - a.se) / (b.fi - a.fi);
}
void Ins(int id, pll pt) {
    int sz = conv[id].size();
    while(sz > 1) {
        pll t1 = conv[id][sz - 2], t2 = conv[id][sz - 1];
        if(Slope(t1, t2) >= Slope(t2, pt)) conv[id].pop_back(), sz--;
        else break;
    } conv[id].pb(pt);
}
ll Getb(int id, ll k) {
    int sz = conv[id].size();
    if(!sz) return inf;
    while(sz > 1) {
        pll t1 = conv[id][0], t2 = conv[id][1];
        if(Slope(t1, t2) <= k) conv[id].pop_front(), sz--;
        else break;
    }
    pll hd = conv[id][0];
    return hd.se - hd.fi * k;
}

int main() {
    n = read(), m = read(), A = read(), B = read(), C = read();
    for(int i = 1; i <= m; i++) {
        x[i] = read(), y[i] = read(), p[i] = read(), q[i] = read();
        tr[(i << 1) - 1] = {i, 1, x[i], p[i]};
        tr[i << 1] = {i, 2, y[i], q[i]};
    } sort(tr + 1, tr + m * 2 + 1);
    Ins(1, {0, 0});
    for(int i = 1; i <= m << 1; i++) {
        Event e = tr[i];
        int id = e.id, st = e.sta, t = e.time;
        if(e.type == 1) {
            ll cost = Getb(st, t);
            if(cost >= inf) {
                res[id] = inf;
                continue;
            }
            cost += C + A * t * t + B * t + t;
            res[id] = cost;
        }
        else {
            ll cost = res[id] + t - p[id];
            if(cost >= inf) continue;
            if(st == n) ans = min(ans, cost);
            Ins(st, {2 * A * t, cost + A * t * t - B * t - t});
        }
    } cout << ans << endl;
    return 0;
}

XIV. P4027 [NOI2007] 货币兑换

设在第 i 天用 d 元能买到 cR_i 数量的 A 劵和 c 数量的 B 劵，有 cR_iA_i+cB_i=d，解得 c=\dfrac d{R_iA_i+B_i}。于是设第 i 天最多能获得多少钱，有：

稍微化简可得 f_i=\max_{1\leq j<i} cR_jA_i+cB_i，其中 c 只与 j 有关。注意到有两个和 i,j 同时有关的项，所以用不起来斜率优化了吗？Nope！由于没有只与 j 有关的项，所以我们可以进行一些移项：将其看作直线 ax+by=c 化简得到 y=\dfrac{c-ax}b，本题中即 c=\dfrac{f_i}{B_i}-cR_j\times \dfrac{A_i}{B_i}。然后使用斜率优化即可。

注意到斜率和插入点横坐标都不单调，所以需要使用李超线段树。虽然查询位置不是整数，但是注意到我们已经知道了所有查询位置 \dfrac{A_i}{B_i}，故将 \dfrac{A_i}{B_i} 离散化即可。

时间复杂度 \mathcal{O}(n\log n)。

const int N = 1e5 + 5;

double d[N], k[N], b[N];
double Get(double x, int id) {return k[id] * x + b[id];}

int mx[N << 2];
void modify(int l, int r, int x, int v) {
    if(l == r) {
        if(Get(d[l], v) > Get(d[l], mx[x])) mx[x] = v;
        return;
    } int m = l + r >> 1;
    if(Get(d[m], v) > Get(d[m], mx[x])) swap(mx[x], v);
    if(Get(d[l], v) > Get(d[l], mx[x])) modify(l, m, x << 1, v);
    else if(Get(d[r], v) > Get(d[r], mx[x])) modify(m + 1, r, x << 1 | 1, v);
}
double query(int l, int r, int p, int x) {
    if(l == r) return Get(d[p], mx[x]);
    int m = l + r >> 1; double ans = Get(d[p], mx[x]);
    if(p <= m) return max(query(l, m, p, x << 1), ans);
    return max(query(m + 1, r, p, x << 1 | 1), ans);
}

int n;
double f[N], A[N], B[N], c[N], R[N];
int main() {
    cin >> n >> f[1];
    for(int i = 1; i <= n; i++)
        cin >> A[i] >> B[i] >> R[i], d[i] = c[i] = A[i] / B[i];
    sort(d + 1, d + n + 1);
    for(int i = 2; i <= n; i++) {
        int j = i - 1; f[i] = f[j];
        double coef = f[j] / (R[j] * A[j] + B[j]);
        k[j] = coef * R[j], b[j] = coef, modify(1, n, 1, j);
        c[i] = lower_bound(d + 1, d + n + 1, c[i]) - d;
        f[i] = max(f[i], query(1, n, c[i], 1) * B[i]);
    } printf("%.3lf\n", f[n]);
    return 0;
}

XV. CF643C Levels and Regions

根据 (1-p)+(1-p)^2+(1-p)^3+\cdots=\dfrac{1}{1-(1-p)}=\dfrac{1}p 可知如果 A 在 B 发生时有 p 的概率发生，那么发生 A 所需的 B 的发生次数的期望为 \dfrac 1p。因此可以设计 DP：

记 s 为 t 的前缀和，方程式改写为：

将 s_k 和 s_j 分离，再记 \dfrac 1 t 的前缀和为 sr，\dfrac s t 的前缀和为 sd，有

显然的斜率优化，时间复杂度 \mathcal{O}(nk)。可不可以 wqs 二分？

const int N = 2e5 + 5;
const ll inf = 1e18;

int n, k, hd, tl, d[N];
double t[N], s[N], p[N], sd[N], sr[N], x[N], y[N], f[N], g[N];
double slope(int i, int j) {return (y[j] - y[i]) / (x[j] - x[i]);}
int main() {
    cin >> n >> k;
    for(int i = 1; i <= n; i++)
        cin >> t[i], s[i] = s[i - 1] + t[i],
        sd[i] = sd[i - 1] + s[i] / t[i], sr[i] = sr[i - 1] + 1 / t[i];
    for(int i = 1; i <= n; i++) g[i] = inf;
    while(k--) {
        d[hd = tl = 0] = 0;
        for(int i = 1; i <= n; i++) {
            while(hd < tl && slope(d[hd], d[hd + 1]) <= sr[i]) hd++;
            int j = d[hd]; f[i] = g[j] + sd[i] - sd[j] - s[j] * (sr[i] - sr[j]);
            x[i] = s[i], y[i] = g[i] - sd[i] + s[i] * sr[i];
            while(hd < tl && slope(d[tl - 1], d[tl]) >= slope(d[tl], i)) tl--; d[++tl] = i;
        } cpy(g, f, N);
    }
    printf("%.10lf\n", f[n]);
    return 0;
}

XVI. CF1083E The Fair Nut and Rectangles

由于所有矩形不会互相包含，所以我们按照 x 递增排序后 y 递减。考虑设计 DP：设 f_i 表示考虑前 i 个矩形的答案且第 i 个矩形必须选，枚举上一个选择的矩形，有：

显然的斜率优化。

不用 fread 擦着时限过，用了只有时限的 1/5，只能说离谱。

const int N = 1e6 + 5;

ll n, f[N], hd, tl, d[N], ans;
struct Rect {
    ll x, y, w;
    bool operator < (const Rect &v) const {return x < v.x;}
} c[N];

ld slope(int i, int j) {return (ld)(f[j] - f[i]) / (c[j].x - c[i].x);}

int main() {
    cin >> n;
    for(int i = 1; i <= n; i++) c[i].x = read(), c[i].y = read(), c[i].w = read();
    sort(c + 1, c + n + 1);
    for(int i = 1; i <= n; i++) {
        while(hd < tl && slope(d[hd], d[hd + 1]) >= c[i].y) hd++;
        int j = d[hd]; f[i] = f[j] + c[i].y * (c[i].x - c[j].x) - c[i].w;
        while(hd < tl && slope(d[tl - 1], d[tl]) <= slope(d[tl], i)) tl--; d[++tl] = i;
        ans = max(ans, f[i]);
    } cout << ans << endl;
    return 0;
}

XVII. LOJ 2769. 「ROI 2017 Day 1」前往大都会

首先求出最短路图 T，是个 DAG。注意到 T 上属于相同的铁路的一段边，我们有转移 f_i=\max_j f_j+(dis_i-dis_j)^2，其中 j 能到达 i。显然的斜率优化。

但是不同于普通斜率优化的是我们要维护一个斜率单调递增的上凸包，这合理吗？没有关系，单调队列不管用一般可以用单调栈，在队尾及时弹出不符合题意的决策即可。

时间复杂度 \mathcal{O}(m\log m)。

const int N = 1e6 + 5;
const ld eps = 1e-9;

int cnt, hd[N], nxt[N], to[N], val[N];
void add(int u, int v, int w) {
    nxt[++cnt] = hd[u], hd[u] = cnt;
    to[cnt] = v, val[cnt] = w;
}

int n, m, v[N], t[N];
struct Edge {
    int u, v, w, id;
};
vector <Edge> e[N];

int dis[N], vis[N], deg[N];
void Dijkstra() {
    priority_queue <pii, vector <pii>, greater <pii>> q;
    mem(dis, 63, N), dis[1] = 0, q.push({0, 1});
    while(!q.empty()) {
        pii t = q.top(); q.pop();
        int id = t.se, ds = t.fi;
        if(vis[id]) continue;
        if(id == n) break;
        vis[id] = 1;
        for(int i = hd[id]; i; i = nxt[i]) {
            int it = to[i], nd = ds + val[i];
            if(dis[it] > nd) q.push({dis[it] = nd, it});
        }
    }
}

vector <int> st[N << 1], v2[N << 1], bel[N];
ll ans[N], r;
ll Y(int x) {return ans[x] + 1ll * dis[x] * dis[x];}
ld slope(int i, int j) {return (ld)(Y(j) - Y(i)) / (dis[j] - dis[i]);}
void checksl(int id, ll sl) {
    int sz = st[id].size();
    while(sz > 1 && slope(st[id][sz - 2], st[id][sz - 1]) <= sl)
        sz--, st[id].pop_back();
}
void checkpt(int id, int p) {
    int sz = st[id].size();
    while(sz > 1 && slope(st[id][sz - 2], st[id][sz - 1]) + eps <= 
        slope(st[id][sz - 1], p)) sz--, st[id].pop_back();
    st[id].pb(p);
}
void dfs(int id) {
    deg[id] = -1;
    map <int, bool> mp; 
    for(int it : bel[id]) {
        while(st[it].size() > 1 && slope(st[it][st[it].size() - 2], st[it].back())
            <= dis[id] * 2) st[it].pop_back();
        if(!st[it].empty()) {
            int tp = st[it].back();
            ans[id] = max(ans[id], ans[tp] + 1ll * (dis[id] - dis[tp]) * (dis[id] - dis[tp]));
        }
    }
    for(int it : bel[id]) checkpt(it, id);
    for(int i = hd[id]; i; i = nxt[i])
        if(dis[id] + val[i] == dis[to[i]] && !--deg[to[i]]) dfs(to[i]);
}
int main() {
    cin >> n >> m;
    for(int i = 1; i <= m; i++) {
        int s = read(); v[1] = read();
        e[i].resize(s);
        for(int j = 1; j <= s; j++)
            t[j] = read(), v[j + 1] = read();
        for(int j = 1; j <= s; j++) {
            add(v[j], v[j + 1], t[j]);
            e[i][j - 1] = {v[j], v[j + 1], t[j], cnt};
        }
    } Dijkstra();
    for(int i = 1; i <= m; i++) {
        v2[++r].pb(e[i][0].u);
        for(Edge it : e[i]) {
            if(dis[it.u] + it.w == dis[it.v]) deg[it.v]++;
            else r++;
            v2[r].pb(it.v);
        }
    }
    for(int i = 1; i <= r; i++)
        for(int it : v2[i]) bel[it].pb(i);
    for(int i = 1; i <= n; i++) if(!deg[i]) dfs(i);
    cout << dis[n] << " " << ans[n] << endl;
    return 0;
}

XVIII. P6173 [USACO16FEB]Circular Barn P

首先看到环不难想到破环成链，枚举断点后变成序列上的问题，DP 状态非常明显：设 f_{i,j} 表示在位置 i 及以前解锁 j 个谷仓且位置 i 下一个位置打开谷仓（即走到 i 的牛不会再往下一个位置走，保证了无后效性）的最小代价，转移为：

把括号拆开，加一个经典前缀和优化，即设 s_i=\sum_{j=1}^i r_i，sd_i=\sum_{j=1}^ir_i\times i，则方程可写为：

斜率优化即可，时间复杂度 \mathcal{O}(n^2k)。感觉可以 wqs 二分（大雾。

const int N=2e3+5;
const ld eps=1e-8;

ll n,k,ans=1e12,r[N],a[N],s[N],sd[N];
ll hd,tl,d[N],x[N],y[N],f[N],g[N];
void solve(){
    x[0]=1,mem(f,0x3f,N),f[0]=0;
    for(int i=1;i<=n;i++)s[i]=s[i-1]+a[i],sd[i]=sd[i-1]+a[i]*i;
    for(int c=1;c<=k;c++){
        hd=tl=1;
        for(int i=1;i<=n;i++){
            while(hd<tl&&y[d[hd+1]]-y[d[hd]]<=(double)s[i]*(x[d[hd+1]]-x[d[hd]]))hd++;
            int j=d[hd]; g[i]=f[j]+sd[i]-sd[j]-(s[i]-s[j])*(j+1);
            x[i]=i+1,y[i]=f[i]-sd[i]+s[i]*(i+1);
            while(hd<tl&&(double)(y[d[tl]]-y[d[tl-1]])*(x[i]-x[d[tl]])>=(double)(y[i]-y[d[tl]])*(x[d[tl]]-x[d[tl-1]]))tl--;
            d[++tl]=i;
        } cpy(f,g,N);
    } ans=min(ans,f[n]);
}
int main(){
    cin>>n>>k;
    for(int i=1;i<=n;i++)cin>>r[i],r[i+n]=r[i];
    for(int i=1;i<=n;i++){
        for(int j=1;j<=n;j++)a[j]=r[i+j-1];
        solve();
    } cout<<ans<<endl;
    return 0;
}

8. 二分队列

8.1 算法介绍

二分队列常用来优化具有决策单调性的 DP 问题，且这里的决策单调性指的是对于 i 和 j 的最优决策点 p_i,p_j，若 i<j 则一定满足 p_i\leq p_j。要求转移时的贡献能快速计算，一般是 \mathcal{O}(1)。

通常情况下的限制为：贡献函数二阶导恒为非负，求最小值 或 二阶导恒为非正，求最大值。

具体地，建立一个存储三元组 (j,l,r) 的队列，表示 l\sim r 的最优决策点是 j。每次判断队首三元组是否过时，若 r<i 则弹出队首。否则将 l 赋值为 i。

加入决策时，如果 i 相比队尾的 j 转移到 l 更优，那么根据决策单调性，i 转移到 l\sim n 比 j 更优，因此 (j,l,r) 就完全无用了，弹出。重复操作直到不满足条件或队列仅剩下一个三元组。

接下来取出队尾的三元组 (j,l,r)，我们要找到一个位置 p 使得 p 以前的位置，从 j 转移更优；而 p 及 p 以后的位置，从 i 转移更优。而这个因为能快速计算贡献，所以显然可以二分出来。注意二分位置从 l 到 r+1。如果 p\leq n 那么将 (i,p,n) 压入队列即可，同时不要忘记将队尾三元组的 r 改为 p-1。

综上，整个算法时间复杂度为 \mathcal{O}(n\log n)。

8.2. 例题

I. P1912 [NOI2009] 诗人小G

经典题。

转移方程为 f_i=\min_{0\leq j<i}|s_i-s_j+(i-j-1)-L|^P，s_i 是长度前缀和。可以使用二分队列。

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

#define ll long long
#define ld long double
#define pb push_back
#define fi first
#define se second
#define mcpy(x,y) memcpy(x,y,sizeof(y))
#define mem(x,v) memset(x,v,sizeof(x))

#define gc getchar()
inline int read(){
    int x=0; char s=gc;
    while(!isdigit(s))s=gc;
    while(isdigit(s))x=x*10+s-'0',s=gc;
    return x;
}

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

const int N=1e5+5;
const int S=30+5;
const ld eps=1e-10;

int n,p,hd,tl,tr[N],s[N];
char str[N][S];
ld f[N],l;

ld cal(int i,int j){return ksm(abs(s[j]-s[i]-l),p);}
ld val(int i,int j){return f[i]+cal(i,j);}

struct tuple{
    int j,l,r;
}d[N];
void print(int x){
    if(!x)return;
    print(tr[x]);
    for(int i=tr[x]+1;i<=x;i++)cout<<(str[i]+1)<<(i==x?'\n':' ');
}

void solve(){
    cin>>n>>l>>p,l++;
    for(int i=1;i<=n;i++){
        scanf("%s",str[i]+1);
        s[i]=s[i-1]+strlen(str[i]+1)+1;
    } d[hd=tl=1]={0,1,n};
    for(int i=1;i<=n;i++){
        while(hd<tl&&d[hd].r<i)hd++; d[hd].l=i;
        int j=d[hd].j; tr[i]=j,f[i]=f[j]+cal(j,i);
        while(hd<tl&&val(i,d[tl].l)<=val(d[tl].j,d[tl].l))tl--;
        int l=d[tl].l,r=d[tl].r+1;
        while(l<r){
            int mid=l+r>>1;
            if(val(i,mid)<=val(d[tl].j,mid))r=mid;
            else l=mid+1;
        } if(l<=n)d[tl].r=l-1,d[++tl]={i,l,n};
    } if(f[n]-eps>1e18)puts("Too hard to arrange");
    else cout<<(ll)(f[n])<<endl,print(n);
    puts("--------------------");
}
int main(){
    int t; cin>>t;
    while(t--)solve();
    return 0;
}

II. P3515 [POI2011]Lightning Conductor

题目中的柿子变一下就是 p=(\max_{i\neq j}a_j+\sqrt{|i-j|})-a_i。正反做一遍二分队列即可。

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

typedef double db;
typedef long long ll;
typedef long double ld;
typedef unsigned long long ull;

#define gc getchar()
#define pb push_back
#define mem(x,v,n) memset(x,v,sizeof(int)*n)
#define cpy(x,y,n) memcpy(x,y,sizeof(int)*n)

const ld Pi=acos(-1);
const ld eps=1e-10;
const ll mod=998244353;

inline int read(){
    int x=0; char s=gc;
    while(!isdigit(s))s=gc;
    while(isdigit(s))x=x*10+s-'0',s=gc;
    return x;
}

const int N=5e5+5;

struct tuple{
    int l,r,k;
}d[N];

int n,hd,tl,a[N];
double f[N],g[N],sqr[N];
double val(int l,int r){return a[l]+sqr[r-l];}
void solve(){
    d[hd=tl=1]={1,n,1},mem(f,0,N*2),f[1]=a[1];
    for(int i=2;i<=n;i++){
        while(hd<tl&&d[hd].r<i)hd++; d[hd].l=i;
        int j=d[hd].k; f[i]=a[j]+sqr[i-j];
        while(hd<tl&&val(i,d[tl].l)-eps>=val(d[tl].k,d[tl].l))tl--;
        int l=d[tl].l,r=d[tl].r+1;
        while(l<r){
            int m=l+r>>1;
            if(val(i,m)-eps>=val(d[tl].k,m))r=m;
            else l=m+1;
        } if(l<=n)d[tl].r=l-1,d[++tl]={l,n,i};
    }
}
int main(){
    n=read();
    for(int i=1;i<=n;i++)a[i]=read(),sqr[i]=sqrt(i);
    solve(),reverse(a+1,a+n+1),cpy(g,f,N*2),solve();
    for(int i=1;i<=n;i++){
        g[i]=max(g[i],f[n-i+1]);
        printf("%d\n",max(0,(int)g[i]+(fabs(g[i]-(int)g[i])>eps)-a[n-i+1]));
    }
    return 0;
}

9. 二分栈

9.1. 算法介绍

二分栈算法常用于有如下决策单调性的 DP 问题中：每个决策点 j 只会被比它更前的决策点 i\ (i<j) 反超。记 v_i 为从决策点 i 转移到当前位置的贡献。

通常情况下的限制为：贡献函数二阶导恒为非负，求最大值 或 二阶导恒为非正，求最小值。这一点可以和前面的二分队列进行对比。

我们可以用一个栈维护可能的决策点，栈顶为当前位置的决策点。考虑如何更新决策点：一个自然的想法是如果栈顶劣于次栈顶，就弹出。但这样是错误的：如果存在 i<j<k 满足 k 优于 j 但 劣于 i，那么这个算法会从 k 而不是 i 转移过来。但是我们有补救的机会：如果 i 反超 j 的时间在 j 反超 k 的时间之前，那么在 i 反超 j 时 j 也不会反超 k 成为最优决策，也就是说决策点 j 已经完全无用了。因此，我们可以在加入决策 k 之前二分出 j 反超 k 的时间 t_1 和 i 反超 j 的时间 t_2，如果 t_1\leq t_2，那么弹出 j。重复上述操作直到栈内只剩 1 个元素或 t_1>t_2，然后压入 k。

不要忘了在每次转移前不断二分出次栈顶反超栈顶的时间 t，如果不大于当前时间 i，那么弹出栈顶。

9.2. 例题

I. P5504 [JSOI2011] 柠檬

经典好题。

设 f_i 为取前 i 只贝壳所能获得的最多柠檬数。一个关键的结论是 DP 仅会在相同颜色的贝壳之间转移，证明不难但不易想到。由于 (st^2)''=2s>0 且求最大值，所以使用二分栈即可。

const int N=1e5+5;

int n,h[N],s[N],p[N];
ll f[N],buc[N];
vector <int> st[N];
ll cal(int p,ll l){return f[p-1]+l*l*s[p];}

#define tp st[c][st[c].size()-1]
#define se st[c][st[c].size()-2]

int chk(int i,int j){
    int l=p[j],r=buc[s[j]]+1;
    while(l<r){
        int m=l+r>>1;
        if(cal(i,m-p[i]+1)<cal(j,m-p[j]+1))l=m+1;
        else r=m;
    } return l;
}

int main(){
    cin>>n;
    for(int i=1;i<=n;i++)s[i]=read(),p[i]=++buc[s[i]];
    for(int i=1;i<=n;i++){
        int c=s[i];
        while(st[c].size()>1&&chk(se,tp)<=chk(tp,i))st[c].pop_back();
        st[c].push_back(i);
        while(st[c].size()>1&&chk(se,tp)<=p[i])st[c].pop_back();
        f[i]=cal(tp,p[i]-p[tp]+1);
    } cout<<f[n]<<endl;
    return 0;
}

10. 整体 DP

10.1. 算法简介

整体 DP 就是用线段树合并维护 DP。

如果看到请催我补一下这个部分的 blog。

10.2. 例题

*I. P4577 [FJOI2018]领导集团问题

设 f_{i,j} 表示以 i 为根的子树 \min w\geq j 的答案，分两种情况讨论：

• 令 f_{i,j}\ (j\leq w_i)\gets\max(f_{i,j},f_{i,w_i}+1)，在合并完儿子后进行该操作。

对于上述转移方程使用整体 DP 即可。

注意到区间 checkmax 的标记永久化无法合并：对于两棵线段树的同一个位置，它的新值应该是两棵树从根到路径的所有节点的标记的 \max 的和而不是所有节点的标记的和的 \max：\max(a_i)+\max(b_i)\neq \max(a_i+b_i)。

由于我们只进行 f_{i,w_i}+1 且 f_{i,j} 是单调递减的，所以注意到每次 checkmax 就相当于区间 +1，标记永久化即可，时间复杂度 \mathcal{O}(n\log n)。

const int N=2e5+5;

int n,node,R[N],ls[N<<5],rs[N<<5],laz[N<<5],mn[N<<5];
void push(int x){mn[x]=min(mn[ls[x]],mn[rs[x]])+laz[x];}
void modify(int l,int r,int ql,int qr,int &x){
    if(!x)x=++node;
    if(ql<=l&&r<=qr)return laz[x]++,mn[x]++,void();
    int m=l+r>>1;
    if(ql<=m)modify(l,m,ql,qr,ls[x]);
    if(m<qr)modify(m+1,r,ql,qr,rs[x]);
    push(x); 
}
int qval(int l,int r,int p,int x){
    if(l==r||!x)return laz[x];
    int m=l+r>>1;
    if(p<=m)return laz[x]+qval(l,m,p,ls[x]);
    return laz[x]+qval(m+1,r,p,rs[x]); 
}
int qpos(int l,int r,int x,int val){
    if(l==r)return l;
    int m=l+r>>1; val-=laz[x];
    if(mn[ls[x]]<=val)return qpos(l,m,ls[x],val);
    return qpos(m+1,r,rs[x],val);
}
int merge(int x,int y){
    if(!x||!y)return x|y;
    ls[x]=merge(ls[x],ls[y]);
    rs[x]=merge(rs[x],rs[y]);
    laz[x]+=laz[y],push(x);
    return x;
}

int w[N],W[N],ans;
vector <int> e[N];
void dfs(int id){
    for(int it:e[id])dfs(it),R[id]=merge(R[id],R[it]);
    int tmp=qval(1,n,w[id],R[id]);
    modify(1,n,qpos(1,n,R[id],tmp),w[id],R[id]);
}
int main(){
    cin>>n;
    for(int i=1;i<=n;i++)cin>>w[i],W[i]=w[i]; sort(W+1,W+n+1);
    for(int i=1;i<=n;i++)w[i]=lower_bound(W+1,W+n+1,w[i])-W;
    for(int i=2,f;i<=n;i++)cin>>f,e[f].pb(i);
    dfs(1),cout<<qval(1,n,1,R[1])<<endl;
    return 0;
}

*II. P6773 [NOI2020] 命运

神仙题。

注意到对于一端在 x 的子树内，另一端在 x 的祖先的限制 (u,v)，若 u 最深的限制被满足，那么所有限制都被满足。因此，设 f_{i,j} 表示一端在以 i 的子树内，另一端在 i 的子树外且 没有被满足的限制 (u,v)的 v 的最大深度为 j 的方案数。

• 若 (x,y) 不选，则有 f_{x,j}\gets\sum_{\max(k,l)=j}f_{x,k}f_{y,l}，即 \sum_{k=0}^j f_{x,j}f_{y,k}+\sum_{k=0}^{j-1}f_{x,k}f_{y,j}
• 若 (x,y) 选，则有 f_{x,j}\gets\sum_{k=0}^{dep_x}f_{x,j}f_{y,k}。

使用前缀和优化后发现子树合并相当于合并两个带乘法标记的线段树。时间复杂度 \mathcal{O}(n\log n)。

下推标记时不需要新建节点，因为没有子节点意味着子节点所表示区间处 DP 值为 0，乘以任何数后仍为 0。

const ll mod = 998244353;
const int N = 5e5 + 5;

void add(ll &x, ll y) {x += y; if(x >= mod) x -= mod;}
void mul(ll &x, ll y) {x = x * y % mod;}

struct Edge {
    int cnt, hd[N], nxt[N << 1], to[N << 1];
    void add(int u, int v) {nxt[++cnt] = hd[u], hd[u] = cnt, to[cnt] = v;}
} edg, lim;

int node, R[N], ls[N << 5], rs[N << 5];
ll sum[N << 5], laz[N << 5];
void pushdown(int x) {
    if(laz[x] != 1) {
        mul(sum[ls[x]], laz[x]), mul(sum[rs[x]], laz[x]);
        mul(laz[ls[x]], laz[x]), mul(laz[rs[x]], laz[x]);
    } laz[x] = 1;
}
void pushup(int x) {sum[x] = (sum[ls[x]] + sum[rs[x]]) % mod;}
void init(int l, int r, int p, int &x) {
    x = ++node, laz[x] = sum[x] = 1;
    if(l == r) return;
    int m = l + r >> 1;
    if(p <= m) init(l, m, p, ls[x]);
    else init(m + 1, r, p, rs[x]);
}
int merge(int l, int r, int x, int y, ll &s1, ll &s2) {
    if(!x && !y) return 0;
    if(!y) return add(s2, sum[x]), mul(laz[x], s1), mul(sum[x], s1), x;
    if(!x) return add(s1, sum[y]), mul(laz[y], s2), mul(sum[y], s2), y;
    if(l == r) {
        ll tmp = sum[x];
        add(s1, sum[y]), mul(sum[x], s1);
        add(sum[x], sum[y] * s2 % mod), add(s2, tmp);
        return x;
    } pushdown(x), pushdown(y);
    int m = l + r >> 1;
    ls[x] = merge(l, m, ls[x], ls[y], s1, s2);
    rs[x] = merge(m + 1, r, rs[x], rs[y], s1, s2);
    return pushup(x), x;
}
int query(int l, int r, int p, int x) {
    if(!x || r <= p) return sum[x];
    ll m = l + r >> 1, ans = 0; pushdown(x);
    if(m < p) ans = query(m + 1, r, p, rs[x]);
    return add(ans, query(l, m, p, ls[x])), ans;
}

int n, m, dep[N];
void dfs(int id, int f) {
    int mxd = 0; dep[id] = dep[f] + 1;
    for(int i = lim.hd[id]; i; i = lim.nxt[i]) mxd = max(mxd, dep[lim.to[i]]);
    init(0, n, mxd, R[id]);
    for(int i = edg.hd[id], it; i; i = edg.nxt[i])
        if((it = edg.to[i]) != f) {
        dfs(it, id);
        ll G = query(0, n, dep[id], R[it]), GG = 0;
        R[id] = merge(0, n, R[id], R[it], G, GG);
    }
}

int main() {
    cin >> n;
    for(int i = 1; i < n; i++) {
        int x = read(), y = read();
        edg.add(x, y), edg.add(y, x);
    } cin >> m;
    for(int i = 1; i <= m; i++) {
        int u = read(), v = read();
        lim.add(v, u);
    } dfs(1, 0), printf("%d\n", query(0, n, 0, R[1]));
    return 0;
}

III. CF490F Treeland Tour

整体 DP 的经典题。对于每个区间维护 LIS 和 LDS，线段树合并直接取 \max 即可。

具体地，设 f/g_{i,j} 分别表示最后一个元素为 j 的 LIS/LDS 最长长度。合并的时候也要更新答案。

时间复杂度 \mathcal{O}(n\log n)。

const int N = 6e3 + 5;
const int L = 1e6;
const int mod = 998244353;

int node, R[N], ls[N << 5], rs[N << 5];
int pre[N << 5], suf[N << 5], ans;
void modify(int l, int r, int p, int v, int &x, int *val) {
    if(!x) x = ++node; val[x] = max(val[x], v);
    if(l == r) return;
    int m = l + r >> 1;
    if(p <= m) modify(l, m, p, v, ls[x], val);
    else modify(m + 1, r, p, v, rs[x], val);
}
int merge(int x, int y) {
    if(!x || !y) return x | y;
    pre[x] = max(pre[x], pre[y]), suf[x] = max(suf[x], suf[y]);
    ans = max(ans, max(pre[ls[x]] + suf[rs[y]], suf[rs[x]] + pre[ls[y]]));
    ls[x] = merge(ls[x], ls[y]), rs[x] = merge(rs[x], rs[y]);
    return x;
}
int query(int l, int r, int ql, int qr, int x, int *val) {
    if(ql > qr) return 0;
    if(!x || ql <= l && r <= qr) return val[x];
    int m = l + r >> 1, ans = 0;
    if(ql <= m) ans = query(l, m, ql, qr, ls[x], val);
    if(m < qr) ans = max(ans, query(m + 1, r, ql, qr, rs[x], val));
    return ans;
}

int n, a[N];
vector <int> e[N];
void dfs(int id, int f) {
    int np = 0, ns = 0;
    for(int it : e[id]) if(it != f) {
        int tpre = 0, tsuf = 0; dfs(it, id);
        tpre = query(1, L, 1, a[id] - 1, R[it], pre);
        tsuf = query(1, L, a[id] + 1, L, R[it], suf);
        ans = max(ans, max(np + tsuf, ns + tpre) + 1);
        np = max(np, tpre), ns = max(ns, tsuf);
        R[id] = merge(R[id], R[it]);
    }
    modify(1, L, a[id], np + 1, R[id], pre);
    modify(1, L, a[id], ns + 1, R[id], suf);
}
int main() {
    cin >> n;
    for(int i = 1; i <= n; i++) cin >> a[i];
    for(int i = 1; i < n; i++) {
        int u, v; cin >> u >> v;
        e[u].pb(v), e[v].pb(u);
    } dfs(1, 0), cout << ans << endl;
    return 0;
}

11. 长链剖分优化 DP

见 简单树论 长链剖分部分。

简要地，长链剖分可以优化树上与深度有关的动态规划。“十二省联考2019 希望” 是一道非常好的入门题，快去试试吧！

12. 轮廓线 DP

轮廓线 DP，又称插头 DP，是一类比较困难的动态规划类型，

12.1. 算法介绍

12.1.1. 概括

要想理解插头 DP，离开例题是不行的。这里简要阐述一下插头 DP 的核心思想：仅记录轮廓线上有用的信息。

实际上，插头 DP 本质上是一种优化过的状压 DP。这类 DP 一般在网格图上进行。不同于状压直接转移一整行，插头是一个一个格子的转移，大大减小了复杂度：n(行数)\times 4^n(转移复杂度)\to n^2(网格数)\times 2^n(转移复杂度)。

12.1.2. 例题

一个经典例题是 P1879 [USACO06NOV]Corn Fields G，虽然它可以用忽略无用状态的状压 DP 通过，但是用来理解插头 DP 再好不过了。

考虑哪些格子的状态对转移有影响：仅有当前格左边和上方的格子。但是我们显然不能仅记录这两个各自的状态因为虽然知道了当前格能否种玉米，但是没办法推出下一个格子上方的状态，因为我们没有记录当前格右上方格子的状态。因此，我们不仅要记录 (r-1,c) 和 (r,c-1)，还要记录 (r,i)\ (i< c) 以及 (r-1,i)\ (i\geq c) 一共 m 个格子的 0/1 状态。而转移的复杂度仅仅需要枚举所有 2^m 种情况，然后 \mathcal{O}(1) 推出下一个状态。

综上，总时间复杂度 \mathcal{O}(n^22^n)。

12.1.3. 总结

所以为什么轮廓线 DP 叫做插头 DP 呢？因为我们将当前格 c 左边和上方与 c 连通的这两个格子内部对 c 有影响的状态称为插头。更一般的，每个格子的状态对于与其相邻（一般是右方和下方）的两个格子的连通性的影响称为插头。这里的 “连通性” 是广义的，可以简要理解为对转移产生的影响。

而由于大部分插头 DP 是逐格转移，因此通常情况下我们有 “右插头” 以及 “下插头”。在例题中，(r-1,c) 的状态是下插头，(r,c-1) 的状态是右插头。

12.2. 扩展与应用

众所周知，插头 DP 常见于与连通性有关的动态规划题目中。很多状压 DP 也可以使用插头 DP 进行优化。这里给出一些常用技巧与注意点：

• 有的时候一个插头的状态可能不止有 / 无两种（例如在哈密顿回路问题中，一个插头可以被表示成左括号或右括号）。当插头种数不是 2 的幂时，提取插头的状态较麻烦而且常数大。为了方便，一般用最接近且大于插头种数的 2 的幂作为状压的进制，例如当种数为 3 时使用 4 进制，种数为 5 时使用 8 进制。

  但是这样会枚举到很多无用状态，例如种数为 5 的 8 进制只要任意一位 >5 就是不合法状态。为了忽略掉无用状态节省空间，可以使用哈希表（如果是连通块相关题目有时还需再加上最小表示法）压缩状态，写起来稍微有一点麻烦（链式前向星）。此外注意清空哈希表时 memset 的复杂度。

• 插头的讨论方法：有 / 无 \times 下 / 右插头，四种情况分别讨论即可。

• 轮廓线上记录 m+1 个状态时，从一行转移至另一行需要将状态整体左移一位（想一想，为什么）。

12.3. 例题

I. P5056 【模板】插头dp

为了区分插头的连通性（防止出现过早出现回路的情况），需要用括号表示法。三进制直接 4 进制哈希即可。轮廓线上记录 m+1 个插头：(r,1)\sim (r,c-1) 的下插头状态，(r,c-1) 的右插头状态以及 (r-1,c)\sim (r-1,m) 的下插头状态，这样才能转移。

• 无下插头，无右插头：由于必须铺，所以新建联通分量：左括号下插头和右括号右插头。
• 无下插头，有右插头：拐弯至当前格下插头或直走至当前格右插头，插头种类显然不变。
• 有下插头，无右插头：拐弯至当前格右插头或直走至当前格下插头，插头种类显然不变。
• 有下插头，有右插头：这种情况就比较麻烦了，需要根据插头种类继续分类讨论：首先显然合并两个插头会使它们湮灭，所以先删掉下插头和右插头。
  • 左括号右插头，右括号下插头：相当于合并一个连通分量分量，合法当且仅当没有别的插头且在最后一个合法格子。
  • 左括号右插头，左括号下插头：将下插头匹配的右括号变成左括号。
  • 右括号右插头，右括号下插头：将右插头匹配的左括号变成右括号。
  • 右括号右插头，左括号下插头：啥都不影响。

此外，插头的延伸需要考虑延伸至的格子是否可以放回路，若不可以，显然转移不合法。另外特殊考虑当前格不合法的情况，只需要原封不动地将状态塞回去即可。

时间复杂度 \mathcal{O}(nm\times 3^m)。常数较小因为哈希表自动帮我们忽略了不合法的情况：有值的状态才会被插入哈希表。

const int N = 15;
const int H = 299987;

struct HashTable {
    int cnt, id[H << 1], nxt[H << 1], hd[H];
    ll val[H << 1];
    void clear() {mem(hd, 0, H), cnt = 0;}
    void insert(int st, ll v) {
        int r = st % H;
        for(int i = hd[r]; i; i = nxt[i])
            if(id[i] == st) return val[i] += v, void();
        val[++cnt] = v, nxt[cnt] = hd[r];
        id[cnt] = st, hd[r] = cnt;
    }
} h[2];
ll ans;
int n, m, edi, edj, now, pre = 1;
int bit[N], bas[N], mp[N][N];
char s;

int main() {
    cin >> n >> m;
    for(int i = 1; i <= n; i++)
        for(int j = 1; j <= m; j++)
            cin >> s, mp[i][j] = s == '.';
    for(int i = n; i && !edi; i--)
        for(int j = m; j && !edj; j--)
            if(mp[i][j]) edi = i, edj = j;
    for(int i = 0; i <= m; i++) bit[i] = i << 1, bas[i] = 1 << bit[i];
    h[0].insert(0, 1);
    for(int i = 1; i <= n; i++) {
        for(int k = 1; k <= h[now].cnt; k++) h[now].id[k] <<= 2;
        for(int j = 1; j <= m; j++) {
            swap(now, pre), h[now].clear();
            for(int k = 1; k <= h[pre].cnt; k++) {
                ll st = h[pre].id[k], dp = h[pre].val[k];
                int lft = st >> bit[j - 1] & 3, up = st >> bit[j] & 3;
                if(!mp[i][j]) {
                    if(!lft && !up) h[now].insert(st, dp);
                } else if(!lft && !up) {
                    if(mp[i + 1][j] && mp[i][j + 1])
                        h[now].insert(st | bas[j - 1] | (bas[j] << 1), dp);
                } else if(lft && !up) {
                    if(mp[i + 1][j]) h[now].insert(st, dp);
                    if(mp[i][j + 1]) h[now].insert(st + lft * (bas[j] - bas[j - 1]), dp);
                } else if(!lft && up) {
                    if(mp[i + 1][j]) h[now].insert(st + up * (bas[j - 1] - bas[j]), dp);
                    if(mp[i][j + 1]) h[now].insert(st, dp);
                }
                else if(lft == 1 && up == 1){
                    int cur = 1;
                    for(int p = j + 1; p <= m; p++) {
                        int t = st >> (p << 1) & 3;
                        cur += t == 1 ? 1 : t == 2 ? -1 : 0;
                        if(cur == 0) {
                            h[now].insert(st - bas[j - 1] - bas[j] - bas[p], dp);
                            break;
                        }
                    }
                } else if(lft == 2 && up == 2) {
                    int cur = -1;
                    for(int p = j - 2; ~p; p--) {
                        int t = st >> (p << 1) & 3;
                        cur += t == 1 ? 1 : t == 2 ? -1 : 0;
                        if(cur == 0) {
                            h[now].insert(st - (bas[j - 1] + bas[j] << 1) + bas[p], dp);
                            break;
                        }
                    }
                } else if(lft == 1 && up == 2) {
                    if(i == edi && j == edj) ans += dp;
                } else if(lft == 2 && up == 1)
                    h[now].insert(st - (bas[j - 1] << 1) - bas[j], dp);
            }
        }
    }
    cout << ans << endl;
    return 0;
}

*II. P4262 [Code+#3]白金元首与莫斯科

一道神仙插头 DP：轮廓线上有 m 个格子，m+1 个插头（定义为向右或向下延伸一个格子），需要记录所有格子的下插头和当前格子左侧的右插头。考虑把每一个位置看做障碍后重新做一遍插头 DP，时间复杂度为 n^2m^22^m，无法接受。GG。

（Trick）接下来就是一个非常神仙的操作了。通常我们在求出一些不可减信息去掉单点后的值时，可以维护一个前缀信息和与后缀信息和，然后快速合并（例如拉格朗日插值在取值连续时，维护前缀积和后缀积避免求逆元从而线性插值）。

对于本题，就是正反各做一遍插头 DP，用空间换时间，即存储每个位置所有插头状态的权值。合并也很有讲究：首先，被挖掉的格子四周不能有插头，而且轮廓线其它地方的插头状态对应，即若 (r,c) 有下插头，则 (r+1,c) 有上插头，因此若确定了一个轮廓线的状态，那么另一个轮廓线与其对应的合法状态是唯一的。因此可以 2^m 合并。

综上，时空复杂度 \mathcal{O}(nm2^m)。

const int N=18;
const int mod=1e9+7;
void add(int &x,int y){x=(x+y)%mod;}

int n,m,mp[N][N],rev[1<<N];
int f[N][N][1<<N],g[N][N][1<<N];
void solve(int f[][N][1<<N]){
    f[1][0][0]=1;
    for(int i=1;i<=n;i++){
        for(int j=1;j<=m;j++)
            for(int k=0;k<1<<m+1;k++)
                if(f[i][j-1][k]){
                    int lft=k>>j-1&1,up=k>>j&1,v=f[i][j-1][k];
                    if(!mp[i][j]){
                        if(!lft&&!up)add(f[i][j][k],v);
                        continue;
                    }
                    if(!lft&&!up){
                        if(mp[i+1][j])add(f[i][j][k^(1<<j-1)],v);
                        if(mp[i][j+1])add(f[i][j][k^(1<<j)],v);
                        add(f[i][j][k],v);
                    }
                    if(!lft&&up)add(f[i][j][k^(1<<j)],v);
                    if(lft&&!up)add(f[i][j][k^(1<<j-1)],v);
                }
        if(i<n)for(int k=0;k<1<<m;k++)
            f[i+1][0][k<<1]=f[i][m][k];
    }
}
int main(){
    cin>>n>>m;
    for(int i=0;i<1<<m+1;i++)for(int j=0;j<m+1;j++)
        rev[i]+=(i>>m-j&1)<<j;
    for(int i=1;i<=n;i++)for(int j=1;j<=m;j++)
        cin>>mp[i][j],mp[i][j]=!mp[i][j];
    solve(f);
    for(int i=1;i<=n;i++)reverse(mp[i]+1,mp[i]+m+1);
    reverse(mp+1,mp+n+1),solve(g);
    for(int i=1;i<=n;i++,cout<<endl)
        for(int j=1;j<=m;j++,cout<<' ')
            if(!mp[n-i+1][m-j+1])cout<<"0";
            else{
                int ans=0;
                for(int k=0;k<1<<m+1;k++){
                    int lft=k>>j-1&1,up=k>>j&1;
                    if(lft||up)continue;
                    add(ans,1ll*f[i][j-1][k]*g[n-i+1][m-j][rev[k]]%mod);
                } cout<<ans;
            }
    return 0;
}

13. DP 套 DP

13.1. 算法简介

DP 套 DP 就是以内层 DP 的结果作为状态进行外层 DP。由于一般动态规划的转移可以看成一个自动机，所以也被称为 DP 自动机。可以类比 AC 自动机上 DP，只是将通过 trie 建出的自动机变成了由内层 DP 建出的自动机。

通常情况下内层 DP 的状态不会很多，否则时间复杂度无法承受。而外层 DP 以内层的状态作为一个维度，因此可以将每个状态重新编号并建出自动机，方便外层 DP 通过自动机转移。

这种类型的动态规划目前在 OI 并不常见。

13.2. 例题

*I. P4590 [TJOI2018]游园会

经典例题，感觉挺神仙的。

这种有奇奇怪怪限制的计数题，从动态规划的角度切入通常比较顺手。显然，DP 的状态设计应与兑奖串的长度（这一维大小为 10^3），以及和 NOI 的匹配长度（这一维大小为 3）有关。

注意到限制 LCS 长度比较棘手，我们考虑求 LCS 的经典方法：f_{i,j}=\max(f_{i-1,j-1}+[a_i=b_j],f_{i-1,j},f_{i,j-1})。 不妨设第一维对应长度为 K 的奖章串 t 的匹配位置，第二维对应兑奖串 s 的匹配位置。一个显然但不易想到的结论是：如果我们知道了 f_{i,j-1}\ (1\leq i\leq K)，根据 s_j 是什么，可以推出所有 f_{i,j}。这启发我们将 LCS 的 DP 数组设计到动态规划里面，转移时根据内层 LCS 的 DP 推出外层 DP 的新状态。此外，不难发现 f_{i,j}\leq f_{i+1,j}\leq f_{i,j}+1，因此其差分数组一定由 0 / 1 组成，这保证我们能够将 f 数组作为 DP 的一维，且大小为 2^K。

时间复杂度 \mathcal{O}(n2^KK(即内层 DP 复杂度){\Sigma}^2)。

template <class T> void cmax(T &a, T b) {a = a > b ? a : b;}
template <class T> void cmin(T &a, T b) {a = a < b ? a : b;}

const int N = 1e3 + 5;
const int K = 15;
const int mod = 1e9 + 7;
void add(int &x, int y) {x = (x + y) % mod;}

int n, k, f[2][1 << K][3], ans[K + 5], s[K + 5];
int main() {
    cin >> n >> k;
    for(int i = 1; i <= k; i++) {
        char tmp; cin >> tmp;
        s[i] = (tmp == 'N' ? 0 : tmp == 'O' ? 1 : 2);
    }
    f[0][0][0] = 1;
    for(int i = 0, p = 0, nw = 1; i < n; i++, swap(p, nw), mem(f[nw], 0, 1 << K))
        for(int j = 0; j < 1 << k; j++)
            for(int m = 0; m < 3; m++) if(f[p][j][m])
                for(int c = 0; c < 3; c++) {
                    int nxt = m == c ? m + 1 : (c == 0 ? 1 : 0), nwS = 0;
                    if(nxt == 3) continue;
                    static int g[K + 2]; g[0] = 0;
                    for(int p = 1; p <= k; p++) g[p] = g[p - 1] + ((j >> p - 1) & 1);
                    for(int p = 1, pre = 0; p <= k; p++) {
                        int cur = max(max(pre, g[p]), g[p - 1] + (c == s[p]));
                        nwS += (cur - pre) << p - 1, pre = cur;
                    } add(f[nw][nwS][nxt], f[p][j][m]);
                }
    for(int i = 0; i < 1 << k; i++) {
        int len = __builtin_popcount(i);
        for(int j = 0; j < 3; j++) add(ans[len], f[n & 1][i][j]);
    }
    for(int i = 0; i <= k; i++) cout << ans[i] << endl;
    return 0;
}

14. 其它技巧

下述优化方法本质上只能算小技巧而非算法，所以合并在一起记录了。

14.1. 双栈优化 DP

如果遇到不可减的信息，但线性次信息加法的时间复杂度可以接受时：维护两个 “对底栈” 并记录栈内信息从栈底到栈顶的前缀和。右端点右移则将新信息压入右边的栈，左端点右移则弹出左边栈顶的信息。若左边的栈为空则将右边的栈的栈顶不断压入左边的栈直到右边的栈为空即可。

不难发现一个信息最多分别进入双栈一次，因此时间复杂度为 \mathcal{O}(nw)，其中 w 是合并信息的复杂度，应用见例题 I。

14.2. 例题

I. 2019 五校联考镇海 B. 小 ω 的仙人掌

题意简述：求最短区间 [l,r] 使得 (w_i,v_i)\ (l\leq i\leq r) 做完背包后权值为 w 的代价 \leq v。

使用对底栈即可，时间复杂度 \mathcal{O}(nw)。写分治 T 掉了，咍咍。

const int N = 1e4 + 5;
const int B = N << 1;
const int W = N >> 1;

int n, ans, wlim, vlim, ltop, rtop, a[N], b[N], e[N];
struct Knapsack {
    int b[W];
    void init() {mem(b, 63, W), b[0] = 0;}
    void ins(int w, int v) {
        for(int i = wlim; i >= w; i--)
            if(b[i - w] + v < b[i])
                b[i] = b[i - w] + v;
    }
} c[N], d[N];
bool check() {
    for(int i = 0; i <= wlim; i++)
        if(c[ltop].b[i] + d[rtop].b[wlim - i] <= vlim)
            return 1;
    return 0;
}
void reverse() {
    while(rtop) {
        ltop++, c[ltop] = c[ltop - 1];
        int id = e[rtop--];
        c[ltop].ins(a[id], b[id]);
    }
}
void push(int id) {
    rtop++, d[rtop] = d[rtop - 1];
    d[rtop].ins(a[id], b[id]), e[rtop] = id;
} 
int main() {
    cin >> n >> wlim >> vlim;
    ans = n + 5, c[0].init(), d[0].init();
    for(int i = 1; i <= n; i++) cin >> a[i] >> b[i];
    for(int i = 1, r = 1; i <= n; i++) {
        while(r <= n && !check()) push(r++);
        if(check()) ans = min(ans, r - i);
        if(ltop == 0) reverse(), ltop--;
        else ltop--;
    } cout << ans << endl;
    return 0;
}