斜率优化初步

· · 算法·理论

首先强调一点斜率优化不一定只用于 dp,它还可以用于很多可以看作斜率的式子中,所以说我们遇到了就具体分析。

What is 斜率优化?

什么是斜率优化?我们来看这样一个问题,现在有一个平面直角坐标系,在其第一象限有很多的点,现在给定一条确定斜率 k(为了简化问题,我们令 k>0)的直线 y=kx+b,在保证这条直线与一个点相交的前提下,使得 b 的取值最小。

很显然,我们可以打暴力,把每个点都去计算一个 b 出来,然后比大小即可,如果是多次询问的话,时间复杂度是 \cal O(n^2) 的,考虑优化。

我们先思考什么样的点一定不可以:

对于图中的 A,B,C 三个点中的 B 点是一定不可以的,从我们的几何直观来看,感觉无论哪种直线,都有可能被 AC 给提前截获,但是为什么呢?我们可以具体证明一下。

Proof:

我们令 A,B,C 的坐标分别为 (x_1,y_1),(x_2,y_2),(x_3,y_3),其中 x_1<x_2<x_3,那么它们之间的几何关系就是 k_{AB}=\frac{y_2-y_1}{x_2-x_1}>k_{BC}=\frac{y_3-y_2}{x_3-x_2}

为了证明 B 是完全无用的,换句话说我们就是要证明对于任意 k 都满足:

很显然,这个式子正向证明是困难的,我们考虑反证。

我们先假设 BA 优,那么它的条件就是 b_B<b_A,把这个式子展开就是 y_2-kx_2<y_1-kx_1,稍微整理一下就是 \frac{y_2-y_1}{x_2-x_1}<k,即 k_{AB}<k

同理,我们假设 BA 优,那么它的条件就是 \frac{y_3-y_2}{x_3-x_2}>k,即 k_{BC}>k

如果说 B 是有用的,那么就应该让 B 同时满足上述的两个条件,即 B 要比 A,B 都更优, 那么其条件就是 k_{AB}<k<k_{BC},但是我们最开始已经规定 k_{AB}>k_{BC},矛盾,故 B 一定无用。

因此我们可以基于此来维护可能为答案的点,我们来观察一下性质:

(图中红线为真正维护的)

不难发现,这样的相邻的两个点之间的斜率是递增的,所以说我们可以用一个单调队列来进行维护,但是如果存储当前点和前一个点之间的斜率的话,是可能会有精度问题的,所以说我们考虑维护每个点的坐标,然后计算斜率时直接用坐标进行交叉相乘然后比较大小即可。

然后考虑对于新加入的一个点,也是直接和单调队列的末端两个点进行上述的比较,如果说末尾的元素满足上述 B 点的情况,就将其弹出,反复这个操作,直到不可为止。最后将新点加入即可。

:::info[Code]{open}

struct Node{ll x,y;}q[INF];
int head=1,tail;

namespace Slope{
    inline void update(Node c){
        while (tail-head>=1){
            Node a=q[tail-1],b=q[tail];
            if ((c.y-b.y)*(b.x-a.x)<=(b.y-a.y)*(c.x-b.x))tail--;
            else break;
        }
        q[++tail]=c;
    }
}

:::

现在的问题就转换为如何快速求的最优点,这个其实也是简单的,我们思考一下最优点的性质。我们假设最优点为 P_i,那么对于这个点是不是一定满足 k_{P_{i-1}P_i}<k<k_{P_iP_{i+1}} 这个条件的?说人话就是最优点和其前一个点所构成的直线的斜率是小于当前直线的斜率的,最优点和其后一个点所构成的直线的斜率是大于当前直线的斜率的。

再加上当前使用的是单调队列,所以说其具有单调性,所以说我们可以采用类似于二分答案的方式进行查找答案。比如说当前枚举到的点是 u,那么就计算 uu-1 的斜率 k_u,然后和当前直线的斜率 k 进行比较,如果说 k_u<k 那么就说明 u 有成为答案的可能性,就先记录下来,然后去 [u+1,r] 这个区间去找,反之就去 [l,u-1] 这个区间去找。

这样时间复杂度就是 \cal O(n\log n) 的了。

:::info[Code]{open}

struct Node{ll x,y;}q[INF];
int head=1,tail;

namespace Slope{
    inline bool check(int a,int b,ll k){return (q[b].y-q[a].y)<=k*(q[b].x-q[a].x);}
    inline int find_best(ll k){
        if (tail==head)return head;
        int l=head+1,r=tail,res=head;
        while (l<=r){
            int mid=(l+r)>>1;
            if (check(mid-1,mid,k))res=mid,l=mid+1;
            else r=mid-1;
        }
        return res;
    }
    inline void update(Node c){
        while (tail-head>=1){
            Node a=q[tail-1],b=q[tail];
            if ((c.y-b.y)*(b.x-a.x)<=(b.y-a.y)*(c.x-b.x))tail--;
            else break;
        }
        q[++tail]=c;
    }
}

:::

当然这里也有一个小小的优化,如果说多条直线的斜率 k 是单调递增的,那么其实就可以不用二分了,可以直接用单调队列进行维护。为什么呢?我们这样考虑:

假设我们有两个决策点 j_1j_2(满足 x_{j_1} < x_{j_2}),且当前查询斜率为 k_i

在我们刚才维护的过程中,如果 j_2j_1 更优,根据之前的推导,必须满足:

K(j_1, j_2) = \frac{y_{j_2} - y_{j_1}}{x_{j_2} - x_{j_1}} \le k_i

因为 k 是单调递增的,所以满足 k_{i+1} \ge k_i

既然 K(j_1, j_2) \le k_i 已经成立,那么显然有:

K(j_1, j_2) \le k_{i+1}

这说明如果 j_2 在当前时刻已经胜过了 j_1,那么在未来斜率更大的时刻,j_2 依然会胜过 j_1

那么 j_1 就一定是无用的了,可以直接弹出,直到不满足上述式子的情况为止。此时的队首就是当前应该取到的点了。

:::info[Code]{open}

struct Node{ll x,y;}q[INF];
int head=1,tail;

namespace Slope{
    inline int get_best(Node c){
        while (tail>head){
            Node a=q[head],b=q[head+1];
            if ((c.y-a.y)*(c.x-b.x)<(c.x-a.x)*(c.y-b.y))head++;
            else break;
        }
        return head;
    }
    inline void update(Node c){
        while (tail-head>=1){
            Node a=q[tail-1],b=q[tail];
            if ((c.y-b.y)*(b.x-a.x)<=(b.y-a.y)*(c.x-b.x))tail--;
            else break;
        }
        q[++tail]=c;
    }
}

::: :::warning{open} 无论是二分还是单调队列,其都满足 y=kx+b 中的 x 单调的性质,如果说 x 不单调,那么就应该考虑使用李超线段树或者 CDQ 分治来做了(斜率优化进阶中会讲)。

横坐标 x 查询斜率 k 维护工具 复杂度
单调 单调 单调队列 (head/tail 双向操作) O(n)
单调 不单调 单调队列 + 二分查找 O(n \log n)
不单调 单调/不单调 李超线段树 / CDQ 分治 O(n \log n)

一个实战技巧:在时间复杂度允许的情况下,其实都可以写李超线段树或者二分,因为其最为完备。 :::

我们回顾上述的点所构成的样子,不难发现这是一个下凸壳,那么有没有可能维护上凸壳呢?这是很显然是有的,维护上凸壳就是使 b 最大,代码中只需要将小于改成大于即可。

这里有初学者经常犯的一个问题,k 的正负对维护上凸壳还是下凸壳有影响吗?这很显然是没有的,你通过画图就可以发现,无论 k 是正还是负,其判断的方式都是一样的,代码是完全没有变化的,所以说与其无关。

最后我们思考一个简单的问题,什么时候维护上凸壳,什么时候维护下凸壳呢?这个答案也是很显然的,通常,我们要求得式子是在 b 这一项,所以说我们根据这个式子是要求最大化还是最小化来判断 b 是应该最大化还是最小化,从而判断是应该维护上凸壳还是下凸壳。

当然这个不是绝对的,有可能我们要求得这个式子就是一个类似于斜率的东西,那么我们就可以选择画出点集和边集,来判断具体是应该维护上凸壳还是下凸壳(下面有例题)。

那么我们回到最初的问题,什么是斜率优化?这个就是斜率优化,通过单调队列在较快的时间内完成对于最小(或最大) b 的求解。那么它如何应用呢?

斜率优化的应用

最常见的应用就是用于优化动态规划,我们以一道例题入手:

例题 1(下凸壳)

P3195 [HNOI2008] 玩具装箱 - 洛谷

我们先思考出暴力的状态转移方程式长成啥样?不难发现,我们设 dp_i 表示以 i 玩具结尾的最小费用,那么有:

dp_i=\min_{0\le j \le i-1} dp_j+[i-(j+1)+\sum_{k=j+1}^{i}C_k-L]^2

其中 i-(j+1)+\sum_{k=j+1}^iC_k 可以用前缀和来求解,即 sum_i-sum_j-1,那么我们就可以做到时间复杂度为 \cal O(n^2) 的暴力,应该有 70

:::info[Code]

#include<bits/stdc++.h>
#define ll long long
using namespace std;
const int INF=5e4+10;

ll n,L,sum[INF],dp[INF];
namespace yixing{
    inline void Main(){
        cin>>n>>L;
        for (int i=1;i<=n;i++){
            int c;cin>>c;
            sum[i]=sum[i-1]+c+1;
        }
        memset(dp,0x3f,sizeof(dp));
        dp[0]=0;
        for (int i=1;i<=n;i++){
            for (int j=0;j<i;j++)dp[i]=min(dp[i],dp[j]+(sum[i]-sum[j]-1-L)*(sum[i]-sum[j]-1-L));
        }
        cout<<dp[n];
    }
}
int main(){
    ios::sync_with_stdio(0);
    cin.tie(0),cout.tie(0);
    yixing::Main();
    return 0;
}

:::

然后我们就考虑优化,我们现在的式子长成这样:

dp_i=\min dp_j+(sum_i-sum_j-1-L)^2

考虑将其展开,同时省去 \min,有:

dp_i=dp_j+sum_i^2+sum_j^2+1+L^2-2\times sum_i \times sum_j-2 \times sum_i-2 \times sum_i \times L+2 \times sum_j+2 \times sum_j\times L+2\times L

我们的斜率优化部分是在对着一条直线做操作,所以说我们这里也应该把它化成一个直线的样子,即化成 y=kx+b 这种形式,我们一个一个来填。

不难发现,当 i 值确定的时候,有很多项都是常数,比如说 sum_i 或者像 L^2,2\times L 这种。结合上述斜率优化部分,我们是每一次去找 b 的最小值,这里我们是去找 dp_i 的最小值,这两个是相似的,所以说 dp_i 一定是放在 b 的这个位置。

同时我们又注意到,这些常数项是一定无法放在 y,k,x 的位置的,所以说我们把所有常数项和 dp_i 都放在 b 的位置,那么就有:

dp_j+sum_j^2-2\times sum_i\times sum_j +2\times sum_j+2\times sum_j \times L=dp_i+2\times sum_i-sum_i^2+2\times sum_i\times L-2\times L-1-L^2

我们现在考虑如何再凑出一个 kx 来。

由于在斜率优化部分,每一条直线的 k 也是确定的,所以说这里合并出来的 k 也应该是一个常数。现在左侧是有两种变量,一个是 dp_j,一个是 sum_j,其中 sum_j 的数量明显大于前者,所以说我们考虑合并 sum_j。不难发现可以将 -2\times sum_i\times sum_j,2\times sum_j,2\times sum_j\times L 合并在一起(不能合并 sum_j^2 因为其提取出来的项还是 sum_j,这个是不确定的)。然后把这个合并了的挪到右侧,作为 kx,那么此时左侧剩下的就是 y 了。

现在这个式子就长成这样:

\underbrace{dp_j+sum_j^2}_{\text{这是 } y}=\underbrace{(2\times sum_i-2-2\times L)}_{\text{这是 }k}\underbrace{sum_j}_{\text{这是 }x}+\underbrace{(dp_i+2\times sum_i-sum_i^2+2\times sum_i\times L-2\times L-1-L^2)}_{\text{这是 } b}

现在所有的项都处理出来了,那么就可以直接套斜率优化的模板了。这里就等价于每一次给定的直线为上面这一条,然后就可以计算出来 dp_i 的值,那么新点的坐标也就有了,为 (sum_i,dp_i+sum_i^2),把这个点加入队列即可。

:::info[Code]

#include<bits/stdc++.h>
#define ll long long
using namespace std;
const int INF=5e4+10;

ll n,L,sum[INF],dp[INF];
struct Node{ll a,b;}q[INF];
int head,tail;

namespace Slope{
    inline bool check(int a,int b,ll k){return (q[b].b-q[a].b)<=k*(q[b].a-q[a].a);}
    inline int get_best(ll k){
        if (head==tail)return head;
        int l=head+1,r=tail,res=head;
        while (l<=r){
            int mid=(l+r)>>1;
            if (check(mid-1,mid,k))res=mid,l=mid+1;
            else r=mid-1;
        }
        return res;
    }
    inline void update(Node c){
        while (tail>head){
            Node a=q[tail-1],b=q[tail];
            if ((b.b-a.b)*(c.a-b.a)>=(b.a-a.a)*(c.b-b.b))tail--;
            else break;
        }
        q[++tail]=c;
    }
}
namespace yixing{
    inline void Main(){
        cin>>n>>L;
        for (int i=1;i<=n;i++){
            int c;cin>>c;
            sum[i]=sum[i-1]+c+1;
        }
        q[++tail]={0,0};
        for (int i=1;i<=n;i++){
            int j=Slope::get_best(2*sum[i]-2*L-2);
            ll x=q[j].a,y=q[j].b;
            ll b=y-(2*sum[i]-2*L-2)*x;
            dp[i]=b-(2*sum[i]+2*sum[i]*L-sum[i]*sum[i]-L*L-2*L-1);
            Node ne_pos={sum[i],dp[i]+sum[i]*sum[i]};
            Slope::update(ne_pos);
        }
        cout<<dp[n];
    }
}
int main(){
    ios::sync_with_stdio(0);
    cin.tie(0),cout.tie(0);
    yixing::Main();
    return 0;
}

:::

这个就是一个标准的斜率优化 dp 的问题,其实上来说,列出基础的状态转移方程,其实就已经做完一大半了,剩下的就几乎是模板了。

我们再来看一道题嘛。

例题 2(上凸壳)

P3628 [APIO2010] 特别行动队 - 洛谷

还是老套路,先列出暴力的状态转移方程,考虑设 dp_i 表示以 i 结尾的队伍修正战斗力的最大值,那么根据题目描述,很显然有:

dp_i=\min _{0\le j \le i-1} dp_j+A\times (\sum_{k=j+1}^i x_k)^2+B\times \sum_{k=j+1}^i x_k+C

其中 \sum_{k=j+1}^i x_k 可以用前缀和维护,所以说有:

dp_i=\min _{0\le j \le i-1} dp_j+A\times (sum_i-sum_j)^2+B\times (sum_i-sum_j)+C

这样做的时间复杂度就是 \cal O(n^2) 的,可以拿 50 分。

:::info[Code]

#include<bits/stdc++.h>
#define ll long long
using namespace std;
const int INF=1e6+10;

ll n,a,b,c;
ll x[INF],sum[INF],dp[INF];
namespace yixing{
    inline void Main(){
        cin>>n>>a>>b>>c;
        for (int i=1;i<=n;i++)cin>>x[i],sum[i]=sum[i-1]+x[i];
        memset(dp,-0x3f,sizeof(dp));
        dp[0]=0;
        for (int i=1;i<=n;i++){
            for (int j=0;j<i;j++){
                ll X=sum[i]-sum[j];
                dp[i]=max(dp[i],dp[j]+(a*X*X+b*X+c));
            }
        }
        cout<<dp[n];
    }
}
int main(){
    ios::sync_with_stdio(0);
    cin.tie(0),cout.tie(0);
    yixing::Main();
    return 0;
}

:::

然后考虑斜率优化,将式子展开有:

dp_i=dp_j+Asum_i^2-2\times Asum_isum_j+Asum_j^2+Bsum_i-Bsum_j+C

简单合并一下并且按照 y=kx+b 的形式分配,有:

\underbrace{dp_j+Asum_j^2-Bsum_j}_{\text{这是 }y}=\underbrace{2\times Asum_i}_{\text{这是 }k}\underbrace{sum_j}_{\text{这是 }x}+\underbrace{(dp_i-Asum_i^2-Bsum_i-C)}_{\text{这是 }b}

然后套斜率优化的模板就可以了。

:::info[Code]

#include<bits/stdc++.h>
#define ll long long
using namespace std;
const int INF=1e6+10;
ll n,a,b,c;
ll x[INF],sum[INF],dp[INF];
struct Node{ll a,b;}q[INF];
int head=1,tail;

namespace Slope{
    inline bool check(int a,int b,ll k){return (q[b].b-q[a].b)>=k*(q[b].a-q[a].a);}
    inline int get_best(ll k){
        if (tail==head)return head;
        int l=head+1,r=tail,res=head;
        while (l<=r){
            int mid=(l+r)>>1;
            if (check(mid-1,mid,k))res=mid,l=mid+1;
            else r=mid-1;
        }
        return res;
    }
    inline void update(Node c){
        while (tail>head){
            Node a=q[tail-1],b=q[tail];
            if ((b.b-a.b)*(c.a-b.a)<=(b.a-a.a)*(c.b-b.b))tail--;
            else break;
        }
        q[++tail]=c;
    }
}
namespace yixing{
    inline void Main(){
        cin>>n>>a>>b>>c;
        for (int i=1;i<=n;i++)cin>>x[i],sum[i]=sum[i-1]+x[i];
        q[++tail]={0,0};
        for (int i=1;i<=n;i++){
            ll k=2*a*sum[i];
            int j=Slope::get_best(k);
            ll x=q[j].a,y=q[j].b;
            ll b1=y-k*x;
            dp[i]=b1+(b*sum[i]+c+a*sum[i]*sum[i]);
            Node ne_p={sum[i],dp[i]+a*sum[i]*sum[i]-b*sum[i]};
            Slope::update(ne_p);
        }
        cout<<dp[n];
    }
}
int main(){
    ios::sync_with_stdio(0);
    cin.tie(0),cout.tie(0);
    yixing::Main();
    return 0;
}

:::

例题 3(多层优化)

斜率优化只能将 \cal O(n^2) 优化成 \cal O(n) 或者 \cal O(n\log n) 吗?显然不是,只要是符合斜率优化的,都可以进行优化,我们看下面这个题目。

P4072 [SDOI2016] 征途 - 洛谷

按照常规,先推暴力,发现这个式子有点难弄,所以说我们考虑对方差做一下操作,有

v^2=\frac{( \overline{v}-v_1)^2+(\overline{v} -v_2)^2+\dots +(\overline{v} -v_m)^2}{m}

展开并合并得:

v^2=\frac{v_1^2+v_2^2+v_3^2+\dots +v_m^2}{m}-\frac{(\sum ^m_{i=1}v_i)^2}{m^2}

然后乘上 m^2 有:

v^2\times m^2=m\times (v_1^2+v_2^2+v_3^2+\dots v_m^2)-(\sum_{i=1}^mv_i)^2

不难发现后者是一个常数,那么我们得目的就是要使前者最小化,那么就很简单了,我们设 dp_{i,l} 表示前 i 段路程分成了 l 天来走的最小平方和。那么就有:

dp_{i,l}=\min _{0\le j\le i-1}dp_{j,l-1}+(\sum _{k=j+1}^{i} x_k)^2

其中后者又可以用前缀和来维护,所以说有:

dp_{i,l}=\min_{0\le j\le i-1} dp_{j,l-1}+(sum_i-sum_j)^2

那么,我们先枚举 l,然后枚举 i,最后枚举 j,用 \cal O(n^3) 的时间复杂度完成暴力,差不多有 80 分的样子。然后我们考虑优化。

:::info[Code]

#include<bits/stdc++.h>
#define ll long long
using namespace std;
const int INF=3e3+10;

ll n,m,v[INF];
ll dp[INF][INF],sum[INF];

namespace yixing{
    inline void Main(){
        cin>>n>>m;
        for (int i=1;i<=n;i++)cin>>v[i],sum[i]=sum[i-1]+v[i];
        memset(dp,0x3f,sizeof(dp));
        ll LIM=dp[0][0];
        dp[0][0]=0;
        for (int k=1;k<=m;k++)for (int i=1;i<=n;i++)for (int j=0;j<i;j++)if (dp[j][k-1]!=LIM)dp[i][k]=min(dp[i][k],dp[j][k-1]+(sum[i]-sum[j])*(sum[i]-sum[j]));
        cout<<dp[n][m]*m-sum[n]*sum[n];
    }
}
int main(){
    ios::sync_with_stdio(0);
    cin.tie(0),cout.tie(0);
    yixing::Main();
    return 0;
}

:::

将上面这个式子展开,并改成 y=kx+b 的形式,有:

dp_{j,l-1}-sum_j^2=2sum_isum_j+dp_{i,l}-sum_i^2

然后就套斜率优化的模板即可。

但是这里要注意一个点,dp_{i,l} 用到的点只有 dp_{j,l-1},所以说 l-1 之前的点都不能用,因为我们每一次都要重新开一个队列来存储,这也是多层优化的一个不同点。

:::info[Code]

#include<bits/stdc++.h>
#define ll long long
using namespace std;
const int INF=3e3+10;

ll n,m,v[INF];
ll dp[INF][INF],sum[INF];
struct Node{ll a,b;}q[INF];
int head=1,tail;

namespace Slope{
    inline bool check(int a,int b,ll k){return (q[b].b-q[a].b)<=k*(q[b].a-q[a].a);}
    inline int get_best(ll k){
        if (head==tail)return head;
        int l=head+1,r=tail,res=head;
        while (l<=r){
            int mid=(l+r)>>1;
            if (check(mid-1,mid,k))res=mid,l=mid+1;
            else r=mid-1;
        }
        return res;
    }
    inline void update(Node c){
        while (tail>head){
            Node a=q[tail-1],b=q[tail];
            if ((b.b-a.b)*(c.a-b.a)>=(b.a-a.a)*(c.b-b.b))tail--;
            else break;
        }
        q[++tail]=c;
    }
}
namespace yixing{
    inline void Main(){
        cin>>n>>m;
        for (int i=1;i<=n;i++)cin>>v[i],sum[i]=sum[i-1]+v[i];
        for (int i=1;i<=n;i++)dp[i][1]=sum[i]*sum[i];
        for (int k=2;k<=m;k++){ 
            head=1,tail=0;
            q[++tail]={0,0};
            for (int i=1;i<=n;i++){
                ll k1=2*sum[i];
                int j=Slope::get_best(k1);
                ll x=q[j].a,y=q[j].b;
                ll b=y-k1*x;
                dp[i][k]=b+sum[i]*sum[i];
                Node ne_p={sum[i],sum[i]*sum[i]+dp[i][k-1]};
                Slope::update(ne_p);
            }
        }
        cout<<dp[n][m]*m-sum[n]*sum[n];
    }
}
int main(){
    ios::sync_with_stdio(0);
    cin.tie(0),cout.tie(0);
    yixing::Main();
    return 0;
}

:::

例题 4(有限制的 j

不难发现,前面的题目,我们枚举的 j 的范围永远都是 [0,i-1],那么如果说 j 的范围不是这个,而是 [0,i-k],其中 k 的值不确定呢?或者说,从前面转移过来的点是有限制性的呢?我们来看这个例题。

P5468 [NOI2019] 回家路线 - 洛谷

还是老套路,先做暴力。

注意到只能对着站点或者某班列车做操作,其中对着站点做操作似乎不太好处理,因为无法处理时间的问题。所以说我们对着某班列车做操作。

我们设 dp_i 表示乘坐第 i 班列车的最小烦躁值,那么我们可以枚举前面乘坐的一班车,然后处理烦躁值,有:

dp_i=\min_{q_j\le p_i\& y_j=x_i} dp_j+A\times (p_i-q_j)^2+B\times (p_i-q_j)+C

那么我们枚举 j,并判断是否满足条件即可在 \cal O(n^2) 的时间复杂度内完成,这样可以拿到 70 分。然后就是要考虑如何优化。

:::info[Code]

#include<bits/stdc++.h>
#define ll long long
using namespace std;
const int MINF=1e6+10,TINF=4e4+10,NINF=1e5+10;

ll n,m,A,B,C;
ll dp[MINF];
struct Line{ll u,v,p,q;}l[MINF];
struct Node{ll a,b;};
vector<Node> q[NINF];
vector<int> t[TINF];
int now_t,head[TINF];

namespace Slope{
    inline int get_best(ll k,int id){
        while (q[id].size()-head[id]>=2){
            Node a=q[id][head[id]],b=q[id][head[id]+1];
            if ((b.b-a.b)<=k*(b.a-a.a))head[id]++;
            else break;
        }
        return head[id];
    }
    inline void update(Node c,int id){
        while (q[id].size()-head[id]>=2){
            Node a=q[id][q[id].size()-2],b=q[id][q[id].size()-1];
            if ((b.b-a.b)*(c.a-b.a)>=(c.b-b.b)*(b.a-a.a))q[id].pop_back();
            else break;
        }
        q[id].push_back(c);
    }
}
namespace yixing{
    bool cmp(const Line a1,const Line a2){return a1.p<a2.p;}
    inline void Main(){
        cin>>n>>m>>A>>B>>C;
        for (int i=1;i<=m;i++)cin>>l[i].u>>l[i].v>>l[i].p>>l[i].q;
        sort(l+1,l+1+m,cmp);
        l[0].u=l[0].p=l[0].q=0,l[0].v=1;
        dp[0]=0;
        ll ans=LONG_LONG_MAX;
        t[0].push_back(0);
        for (int i=1;i<=m;i++){
            while (now_t<=l[i].p){
                for (int j:t[now_t]){
                    int id=l[j].v;
                    Node c={l[j].q,dp[j]+A*l[j].q*l[j].q-B*l[j].q};
                    Slope::update(c,id);
                }
                now_t++;
            }
            int u=l[i].u;
            if ((int)q[u].size()>head[u]){
                ll k=2*l[i].p*A;
                int j=Slope::get_best(k,u);
                ll x=q[u][j].a,y=q[u][j].b;
                ll b=y-k*x;
                dp[i]=b+(B*l[i].p+A*l[i].p*l[i].p+C);
                t[l[i].q].push_back(i);
                if (l[i].v==n)ans=min(ans,dp[i]+l[i].q);
            }else dp[i]=LONG_LONG_MAX;

        }
        cout<<ans;
    }
}
int main(){
    ios::sync_with_stdio(0);
    cin.tie(0),cout.tie(0);
    yixing::Main();
    return 0;
}

:::

我们还是按照常规对其进行拆开,有:

dp_i=dp_j+Ap_i^2-2\times Ap_iq_j+Aq_j^2+Bp_i-Bq_j+C

然后稍加整理,有:

dp_j+Aq_j^2-Bq_j=2\times Ap_iq_j+dp_i-Ap_i^2-Bp_i-C

但是这里很显然是不能直接套用斜率优化的,因为 j 具有两个限制条件,我们遇到这种限制条件,首选的方法就是分组。就像例题 3 一样,每一次循环,都会重新队列清空,那么这里也是,能不能直接对着每一班列车开一个队列呢?显然可以,但是这样就无法判断是否满足 y_j=x_i 这个条件了。

那么我们就考虑对着每一个站点开一个队列,用来维护到达这个站点的列车(即 v_i 等于这个站点)的斜率情况。

最后就是要处理时间的问题了,这里也是有一个非常经典的操作,就是延迟入队。什么意思呢,简单解释来说就是计算出来 dp_i 的值的时候,先不要将其入队,等到某个时候再进行入队。

回顾我们上面所写的代码,几乎都是计算完一个 dp_i 之后马上就将其丢到队列里面去,但是这里多了一个限制条件,时间。换句话说,对于第 i 班列车,等到其到站的时间是 q_i,但是对于我们所枚举的第 i+1 班列车,其出发时间 p_{i+1} 可能满足 p_{i+1}<q_i 的情况,即第 i 班列车还没有到达,第 i+1 班列车已经出发。

那么如果说将 i 加入了队列的话,就有可能会使 i+1 调用到 i 的答案,这个很显然是不正确的转移,从而使答案错误。

那么我们就考虑如何限制一下时间:对于第 i 班列车,其可用的点是不是一定满足 t\le p_i,为了不出现删点又加点的情况,所以说我们考虑在最开始就对着 p_i 进行一次排序。

接下来我们就应该考虑如何进行延迟加点的问题了,很显然,我们可以开一个变量 now 用来表示当前的时间,同时开一个动态数组对于每一个时间点的点进行存储(因为 q 相同的点可能会有很多个,同时他们不一定是以一个先后的顺序出现的,所以说不能只开一个动态数组用来维护这些点),最后每到一个新点,我们就将 now 跳到这个位置,同时在跳的过程中将点加入对应的队列中。

注意到每个点还是只会入队出队一次,所以说时间复杂度是没有变化的。

:::info[Code]

#include<bits/stdc++.h>
#define ll long long
using namespace std;
const int INF=1e6+10;

ll n,m,A,B,C;
ll dp[INF];
struct Line{ll u,v,p,q;}l[INF];
struct Node{ll a,b;};
vector<Node> q[INF];
vector<int> t[INF];
int now_t,head[INF];

namespace Slope{
    inline int get_best(ll k,int id){
        while (q[id].size()-head[id]>=2){
            Node a=q[id][head[id]],b=q[id][head[id]+1];
            if ((b.b-a.b)<=k*(b.a-a.a))head[id]++;
            else break;
        }
        return head[id];
    }
    inline void update(Node c,int id){
        while (q[id].size()-head[id]>=2){
            Node a=q[id][q[id].size()-2],b=q[id][q[id].size()-1];
            if ((b.b-a.b)*(c.a-b.a)>=(c.b-b.b)*(b.a-a.a))q[id].pop_back();
            else break;
        }
        q[id].push_back(c);
    }
}
namespace yixing{
    bool cmp(const Line a1,const Line a2){return a1.p<a2.p;}
    inline void Main(){
        cin>>n>>m>>A>>B>>C;
        for (int i=1;i<=m;i++)cin>>l[i].u>>l[i].v>>l[i].p>>l[i].q;
        sort(l+1,l+1+m,cmp);
        l[0].u=l[0].p=l[0].q=0,l[0].v=1;
        dp[0]=0;
        ll ans=LONG_LONG_MAX;
        t[0].push_back(0);
        for (int i=1;i<=m;i++){
            while (now_t<=l[i].p){
                for (int j:t[now_t]){
                    int id=l[j].v;
                    Node c={l[j].q,dp[j]+A*l[j].q*l[j].q-B*l[j].q};
                    Slope::update(c,id);
                }
                now_t++;
            }
            int u=l[i].u;
            if ((int)q[u].size()>head[u]){
                ll k=2*l[i].p*A;
                int j=Slope::get_best(k,u);
                ll x=q[u][j].a,y=q[u][j].b;
                ll b=y-k*x;
                dp[i]=b+(B*l[i].p+A*l[i].p*l[i].p+C);
                t[l[i].q].push_back(i);
                if (l[i].v==n)ans=min(ans,dp[i]+l[i].q);
            }else dp[i]=LONG_LONG_MAX;

        }
        cout<<ans;
    }
}
int main(){
    ios::sync_with_stdio(0);
    cin.tie(0),cout.tie(0);
    yixing::Main();
    return 0;
}

:::

例题 5(斜率式子的转化)

上述的题目全部都是一个套路,推出式子,然后化成 y=kx+b 的形式,然后用斜率优化。那么如果说这个式子无法化成 y=kx+b 的形式呢?我们来看这道题。

P1721 [NOI2016] 国王饮水记 - 洛谷

这道题毕竟还是一道黑题,肯定没有那么简单,所以说我们一步一步来分析。

Part 1:什么会对答案造成影响?

我们先思考什么样的水箱会对答案造成影响?不难发现应该是 h_i>h_1 的水箱才有贡献的可能,而小于等于的水箱所带来的贡献肯定是小于等于 0 的。

那么我们可以在输入的时候把所有 h_i\le h_1 的水箱全部舍去,只留下对答案有贡献的。

Part 2:怎么操作才可以使的答案最大?

假设说现在有高度为 h_1<h_i<h_j 的三个水箱,考虑分类讨论。

如果说先合并 h_1h_i,然后再合并 h_j 的话,最终的高度是 \frac{1}{2}\times (\frac{1}{2}\times (h_1+h_i)+h_j)=\frac 1 4h_1+\frac 1 4h_i+\frac 1 2h_j

如果说先合并 h_1h_j,然后再合并 h_i 的话,最终的高度是 \frac{1}{2}\times (\frac{1}{2}\times (h_1+h_j)+h_i)=\frac 1 4h_1+\frac 1 4h_j+\frac 1 2h_i

很显然,前者大于后者,所以说一定是先合并矮的,再去合并高的,这样才可以使得答案最大化。

所以说我们可以把 h 按照从小到大的顺序进行排序,然后再进行合并操作。

Part 3:怎么利用动态规划解决?

由于说我们现在已经完成了对 h 的操作,那么接下来的动态规划一定也是从小到大进行递推。类似于征途那道题,我们设 dp_{i,l} 表示前 i 个水箱,分成 l 次合并所能使 h_1 到达的最大高度。根据 \overline h^{'}=\frac{\overline{h}+\sum_{k=j+1}^{i}h_k}{i-j+1} 有:

dp_{i,l}=\max_{0\le j<i} \frac{dp_{j,l-1}+\sum_{k=j+1}^{i}h_k}{i-j+1}

其中 \sum_{k=j+1}^{i}h_k 可用前缀和进行替换,有:

dp_{i,l}=\max_{0\le j<i} \frac{dp_{j,l-1}+sum_i-sum_j}{i-j+1}

这样写的时间复杂度是 \cal O(n^2k) 的,大概能有 32 分。

Part 4:怎么优化动态规划?(关键点)

很显然这个式子我们是没有办法像之前一样进行划分的,但是我们回想一下斜率的定义式:k=\frac{y_2-y_1}{x_2-x_1},不难发现这个玩意似乎是一个除法。

那么能不能也这样处理一下呢?显然是可以的。

这个式子的右侧可以分一下组,看成:

dp_{i,l}=\max_{0\le j<i} \frac{sum_i-(sum_j-dp_{j,l-1})}{i-(j-1)}

其中,有关 i 的项全部都是已知数,有关 j 的项是不确定的,所以说考虑看作两个点 M(i,sum_i),N(j-1,sum_j-dp_{j,l-1}),而 dp_{i,l} 就是这两个点的斜率。

那么现在的问题就转化为了,给定一个确定的点 M,从一个点集 \cal P 中找出一个点 N 使得二者的斜率最大。

由于 i 恒大于 j-1sum_i 恒大于 sum_j-dp_{j,l-1},所以说有点 M 恒在 N 的右上方,又要使得二者的斜率最大,不难想到维护下凸壳。

又注意到 isum_i 都是单调递增的,所以说 M 也一定是在往右上方走,所以说可以直接用单调队列直接优化,而不需要二分,省去了一个 \log n,这样写时间复杂度差不多是 \cal O(nk) 的,但是在 k\le 10^9 这个级别还是过不了。

Part 5:怎么优化掉 k

首先有一个较为显然的事实:每一次合并,至少都要包含一个之前未连通的城市,否则水位不会上升。而除了首都之外只剩下 n-1 个城市了,所以说将 kn 二者取一个 \min 即可,现在这样的时间复杂度是 \cal O(n^2) 的,显然是可以通过的了,但是你会发现,代码交上去只有 55 分左右,剩下错的点,大体上都计算正确了,但是还有有几百或几十的误差,为什么呢?

注意到,我们不可能每一次都去使用题目中给的高精度小数来操作(因为操作一次的时间复杂度至少都是 \cal O(P) 的),所以说无论是 dp 还是 sum 我们都只能用 long double 来存,如果经历 \cal O(n^2) 次上述的状态转移,很显然精度是要炸的,那么还要考虑优化。

Part 6:怎么处理精度爆炸的问题?

设初始状态下,水位与最高城市的差距最大为 H = n \cdot \max(h_i)\Delta=\min_i{h_i-h_{i-1}}L 为当前合并区间的长度。

初始距离很显然为 D_{start} \approx O(nh),当距离 D 缩小到小于 \Delta 时,由于 h_i - h_{i-1} \ge \Delta,此时水位 f 已经非常接近 h_i。在这种情况下,任何 L > 1 的合并都会因为“加入了一个水位更低的城市”而导致最终水位不如 L=1 时的贪心结果。

因此,合并次数 m 满足:

( \text{收缩因子} )^m \cdot (nh) \le \Delta

解得:

m = O(\log \frac{nh}{\Delta})

我们考虑一下极限情况,有 n=8000,h=10^5,\Delta=12.5,此时有 m=\log_2 \frac{nh}{\Delta}\approx 25。那么很显然就可以将 k25 取一个 \min,记作 lim。然后每一次转移都是记录从哪里转移过来的,然后统一用一个高精度小数来重新模拟一遍,然后再将剩下的次数给使用殆尽。

然后此时你再去交,你就会发现,可以拿大概 73 分。

为什么呢?还是那句话,精度要炸,尽管这样已经提升的很大一部分的精度问题了,但是请注意,即使是在 Linux 环境下,long double 的小数位数也只有 18\sim 19 位,所以说我们应该将 k18 取一个 \min,从而保证没有精度的问题,最后再来跑。

删去了高精度小数的部分,只展示核心代码。

:::info[Code]

#define ll long long
#define ld long double
const int INF=8100;
using namespace std;
Decimal ans;
int n,k,p,h[INF],tot=1;//tot 表示有多少个数大于 h1
int pos[INF][20];//pos_il 表示在 i 这个位置,第 l 次操作了的点
ll sum[INF];
ld dp[INF][20];
struct Node{ld a,b;}q[INF];
int head=1,tail=0,id[INF];

namespace Slope{
    inline int get_best(Node c){
        while (tail>head){
            Node a=q[id[head]],b=q[id[head+1]];
            if ((c.b-a.b)*(c.a-b.a)<(c.a-a.a)*(c.b-b.b))head++;
            else break;
        }
        return head;
    }
    inline void udpate(Node c,int p){
        while (tail>head){
            Node a=q[id[tail-1]],b=q[id[tail]];
            if ((b.b-a.b)*(c.a-b.a)>(b.a-a.a)*(c.b-b.b))tail--;
            else break;
        }
        id[++tail]=p;
    }
}
namespace yixing{
    Decimal calc(int city,int tim){
        if (!tim)return h[1];
        return ((calc(pos[city][tim],tim-1)+(sum[city]-sum[pos[city][tim]]))/(city-pos[city][tim]+1));
    }
    inline void Main(){
        cin>>n>>k>>p>>h[1];
        for (int i=2;i<=n;i++){
            int tmp;cin>>tmp;
            if (tmp>h[1])h[++tot]=tmp;
        }
        n=tot;
        sort(h+1,h+1+n);
        for (int i=1;i<=n;i++)sum[i]=sum[i-1]+h[i];
        for (int i=1;i<=n;i++)dp[i][0]=h[1];
        k=min(k,n);
        int lim=min(k,18);
        for (int l=1;l<=lim;l++){
            head=1,tail=0;
            id[++tail]=1;
            for (int i=1;i<=n;i++)q[i]={1.0*(i-1),(sum[i]-dp[i][l-1])};
            for (int i=2;i<=n;i++){
                Node u={1.0*i,1.0*sum[i]};
                int j=id[Slope::get_best(u)];
                pos[i][l]=j;
                dp[i][l]=(sum[i]-(sum[j]-dp[j][l-1]))/(i-(j-1));
                Slope::udpate(q[i],i);
            }
        }
        int rest=n-k+lim,tim;
        ld maxn=0;
        for (int l=0;l<=lim;l++)if (dp[rest][l]>maxn)maxn=dp[rest][l],tim=l;
        ans=calc(rest,tim);
        for (int i=rest+1;i<=n;i++)ans=(ans+h[i])/2;
        cout<<ans.to_string(2*p);
    }
}
int main(){
    ios::sync_with_stdio(0);
    cin.tie(0),cout.tie(0);
    yixing::Main();
    return 0;
}

:::