抱歉，您的浏览器无法访问本站
本页面需要浏览器支持（启用）JavaScript
了解详情 >

CDQ分治优化斜率优化

发布于：2019年8月20日

例题1：

P2365 任务安排

我知道你们会n2n^2大暴力，也会O(n)O(n)普通斜率优化，但是不妨提高一下我们的姿势水平，假设1≤N≤1000001 \le N \le 100000，而且TiT_i和CiC_i不一定是正整数。

朴素dpdp方程：dp[i]=dp[j]−(s+sumt[i])×sumc[j]+sumt[i]×sumc[i]+s×sumc[n]dp[i]=dp[j]-(s+sumt[i]) \times sumc[j]+sumt[i]\times sumc[i]+s\times sumc[n]

x=sumc[j],y=dp[j],k=s+sumt[i],b=sumt[i]×sumc[i]+s×sumc[n]x=sumc[j],y=dp[j],k=s+sumt[i],b=sumt[i]\times sumc[i]+s\times sumc[n]

式子化成dp[i]=y−kx+bdp[i]=y-kx+b的形式。

注意到sumt,sumcsumt,sumc不一定单调，于是我们需要动态地在凸壳上面插入和查询点，但是这样写太过naive了，考虑离线，使用CDQ分治计算答案，每次我们建出左半边的凸壳，查询右半边凸壳的答案，伪代码如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
procedure CDQ(l,r){
    if(l==r){
        用a[l]对应的dp值来更新a的(x,y)
        return
    }
    mid=(l+r)/2
    把[l,r]内的元素按id分成两个区间，id<=mid的分到[l,mid],否则分到[mid+1,r]
    CDQ(l,mid),先递归计算id值小的
    把[l,mid]区间对应的点建出凸壳，用于下一步更新[mid+1,r]的dp值
    for(k in [mid+1,r]){
        用单调队列更新a[k]对应的dp值，因为我们cdq之前按斜率排过序，所以不会出问题
        注意到在这里面即使a[k]对应的dp值的最优决策不在[mid+1,r]中，那么回溯到一个更大的区间后，也会被更新到
    }
    CDQ(mid+1,r)，继续递归下去更新
    按x,y归并排序[l,mid],[mid+1,r],这样回溯的时候x单调递增，才可以直接单调栈(队列)建凸壳
}
int main(){
    输入
    按照斜率对点排序
}

摘自巨佬ypy的博客

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
#include <bits/stdc++.h>
#define MAXN 300005
#define int long long
using namespace std;
inline int read(){
	int x=0,f=1;
	char ch=getchar();
	while (ch<'0'||ch>'9'){
		if (ch=='-') f=-1;
		ch=getchar();
	}
	while (ch>='0'&&ch<='9'){
		x=(x<<3)+(x<<1)+(ch^'0');
		ch=getchar();
	}
	return x*f;
}
int N,S;
int sumT[MAXN],sumC[MAXN];
int ans[MAXN];
struct Point{
	double x,y;
	int id;
}p[MAXN],t[MAXN],q[MAXN];
inline bool cmp1(const Point &A,const Point &B){
	return sumT[A.id]<sumT[B.id];
}
inline bool cmp2(const Point &A,const Point &B){
	if (A.x!=B.x) return A.x<B.x;
	else return A.y<B.y;
}
inline double operator * (const Point &A,const Point &B){
	return A.x*B.y-A.y*B.x;
}
inline Point operator - (const Point &A,const Point &B){
	return Point{A.x-B.x,A.y-B.y};
}
inline double Slope(const Point &A,const Point &B){
	return (double)(A.y-B.y)/(double)(A.x-B.x);
}
inline void Msort1(int l,int r){
	int mid=(l+r)>>1;
	int j=l,k=mid+1;
	for (register int i=l;i<=r;++i){
		if (p[i].id<=mid) t[j++]=p[i];
		else t[k++]=p[i];
	}
	for (register int i=l;i<=r;++i) p[i]=t[i];
}
inline void Msort2(int l,int r){
	int mid=(l+r)>>1;
	int j=l,k=mid+1;
	for (register int i=l;i<=r;++i){
		if (j<=mid&&(k>r||cmp2(p[j],p[k]))) t[j++]=p[i];例题
		else t[k++]=p[i];
	}
	for (register int i=l;i<=r;++i) p[i]=t[i];
}
void CDQ(int l,int r){
	if (l==r){
		p[l].x=(double)sumC[p[l].id];
		p[l].y=(double)ans[p[l].id];
		return ;
	}
	int mid=(l+r)>>1;
	Msort1(l,r);
	CDQ(l,mid);
	int head=1,rear=0;
	for (register int i=l;i<=mid;++i){
		while (head<rear&&((q[rear]-q[rear-1])*(p[i]-q[rear-1]))<0) rear--;
		q[++rear]=p[i];
	}
	for (register int k=mid+1;k<=r;++k){
		int i=p[k].id;
		while (head<rear&&Slope(q[head+1],q[head])<(S+sumT[i])) head++;
		int j=q[head].id;
		ans[i]=min(ans[i],ans[j]+sumT[i]*(sumC[i]-sumC[j])+S*(sumC[N]-sumC[j]));
	}
	CDQ(mid+1,r);
	Msort2(l,r);
}
#undef int
int main(){
#define int long long
	N=read(),S=read();
	for (register int i=1;i<=N;++i){
		sumT[i]=sumT[i-1]+read();
		sumC[i]=sumC[i-1]+read();
		p[i].id=i;
	}
	sort(p+1,p+1+N,cmp1);
	memset(ans,0x3f,sizeof(ans));
	ans[0]=0;
	CDQ(0,N);
	printf("%lld\n",ans[N]);
}

例题2：

P4655 [CEOI2017]Building Bridges

很容易得出dp方程：

dp[i]=min(dp[i],dp[j]+(h[i]−h[j])2+w[i−1]−w[j])dp[i]=min(dp[i],dp[j]+(h[i]-h[j])^2+w[i-1]-w[j])其中w[i]w[i]是前缀和。

考虑斜率优化，化成dp[i]=dp[j]+h[j]2−w[j]−2h[i]h[j]+h[i]2dp[i]=dp[j]+h[j]^2-w[j]-2h[i]h[j]+h[i]^2的形式，设y=dp[j]+h[j]2−w[j]y=dp[j]+h[j]^2-w[j]的，x=h[j]x=h[j]，k=2h[i]k=2h[i]，b=h[i]2b=h[i]^2

注意到h[j]h[j]不一定单调，所以需要CDQ分治。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
#include <bits/stdc++.h>
#define MAXN 100005
#define int long long
using namespace std;
inline int read(){
    int x=0,f=1;
    char ch=getchar();
    while (ch<'0'||ch>'9'){
        if (ch=='-') f=-1;
        ch=getchar();
    }
    while (ch>='0'&&ch<='9'){
        x=(x<<3)+(x<<1)+(ch^'0');
        ch=getchar();
    }
    return x*f;
}
struct Point{
    double x,y;
    int i;
}p[MAXN],t[MAXN],q[MAXN];
inline double operator * (const Point &A,const Point &B){
    return A.x*B.y-A.y*B.x;
}
inline Point operator - (const Point &A,const Point &B){
    return Point{A.x-B.x,A.y-B.y};
}
inline bool cmp1(const Point &A,const Point &B){
    if (A.x!=B.x) return A.x<B.x;
    return A.y>B.y;
}
int h[MAXN],sum[MAXN],dp[MAXN];
inline bool cmp2(const Point &A,const Point &B){
    return h[A.i]<h[B.i];
}
inline void Merge1(int l,int r){
    int mid=(l+r)>>1;
    int j=l,k=mid+1;
    for (register int i=l;i<=r;++i){
        if (p[i].i<=mid) t[j++]=p[i];
        else t[k++]=p[i];
    }
    for (register int i=l;i<=r;++i) p[i]=t[i];
}
inline void Merge2(int l,int r){
    int mid=(l+r)>>1;
    int j=l,k=mid+1;
    for (register int i=l;i<=r;++i){
        if (j<=mid&&(k>r||cmp1(p[j],p[k]))) t[i]=p[j++];
        else t[i]=p[k++];
    }
    for (register int i=l;i<=r;++i) p[i]=t[i];
}
void CDQ(int l,int r){
    if (l==r){
        p[l].x=(double)h[l];
        p[l].y=(double)dp[l]-sum[l]+h[l]*h[l];
        return ;
    }
    int mid=(l+r)>>1;
    Merge1(l,r);
    CDQ(l,mid);
    int head=1,rear=0;
    for (register int i=l;i<=mid;++i){
        while (head<rear&&(q[rear]-q[rear-1])*(p[i]-q[rear-1])<=0) rear--;
        q[++rear]=p[i];
    }
    for (register int k=mid+1;k<=r;++k){
        int i=p[k].i;
        while (head<rear&&(q[head+1]-q[head])*Point{1,(double)2*h[i]}>=0) head++;
        int j=q[head].i;
        dp[i]=min(dp[i],dp[j]+(h[i]-h[j])*(h[i]-h[j])+sum[i-1]-sum[j]);
    }
    CDQ(mid+1,r);
    Merge2(l,r);
}
#undef int
int main(){
#define int long long
    int n=read();
    for (register int i=1;i<=n;++i) h[i]=read(),p[i].i=i;
    for (register int i=1;i<=n;++i) sum[i]=sum[i-1]+read();
    memset(dp,0x3f,sizeof(dp));
    dp[1]=0;
    sort(p+1,p+1+n,cmp2);
    CDQ(1,n);
    printf("%lld\n",dp[n]);
}

例题3：

P4027 [NOI2007]货币兑换

注意到题目有这么一句话：

必然存在一种最优的买卖方案满足：

每次买进操作使用完所有的人民币；

每次卖出操作卖出所有的金券。

所以我们列出dp方程，F[i]F[i]表示第ii天最多的钱数

F[i]=max(F[i−1],F[j]/(A[j]×R[j]+B[j])×(A[i]×R[j]+B[i]))F[i]=max(F[i-1],F[j]/(A[j] \times R[j]+B[j])\times (A[i]\times R[j]+B[i]))

F[i]=F[i−1]F[i]=F[i-1]表示第ii天什么都不做，

F[i]=F[j]/(A[j]×R[j]+B[j])×(A[i]×R[j]+B[i])F[i]=F[j]/(A[j] \times R[j]+B[j]) \times (A[i] \times R[j]+B[i])

表示把第jj天的钱全部按照R[j]R[j]的比例买进股票，然后在第ii天按照R[j]R[j]的比例全部卖出。

考虑斜率优化，由于F[j]/(A[j]×R[j]+B[j])F[j]/(A[j] \times R[j]+B[j])里面只有jj，不妨考虑换元。

设F[j]/(A[j]×R[j]+B[j])=TF[j]/(A[j] \times R[j]+B[j])=T，那么我们有F[i]=T×(A[i]×R[j]+B[i])F[i]=T \times (A[i] \times R[j]+B[i])

拆开，得F[i]=T×A[i]×R[j]+T×B[i]F[i]=T \times A[i] \times R[j]+T \times B[i]

令y=T×R[j]y=T \times R[j]，x=Tx=T，k=−B[i]/A[i]k=-B[i]/A[i]

直接套用CDQ分治模板。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
#include <bits/stdc++.h>
#define MAXN 300005
#define int long long
#define INF 1e12
using namespace std;
inline int read(){
	int x=0,f=1;
	char ch=getchar();
	while (ch<'0'||ch>'9'){
		if (ch=='-') f=-1;
		ch=getchar();
	}
	while (ch>='0'&&ch<='9'){
		x=(x<<3)+(x<<1)+(ch^'0');
		ch=getchar();
	}
	return x*f;
}
int N;
double S,A[MAXN],B[MAXN],R[MAXN],ans[MAXN];
struct Point{
	double x,y;
	int id;
}p[MAXN],t[MAXN],q[MAXN];
inline bool cmp1(const Point &X,const Point &Y){
	return -B[X.id]/A[X.id]>-B[Y.id]/A[Y.id];
}
inline bool cmp2(const Point &A,const Point &B){
	if (A.x!=B.x) return A.x<B.x;
	else return A.y>B.y;
}
inline double operator * (const Point &A,const Point &B){
	return A.x*B.y-A.y*B.x;
}
inline Point operator - (const Point &A,const Point &B){
	return {A.x-B.x,A.y-B.y};
}
inline Point operator + (const Point &A,const Point &B){
	return {A.x+B.x,A.y+B.y};
}
inline double Slope(const Point &A){
	if (A.x==0) return A.y>0?INF:-INF;
	return A.y/A.x;
}
inline void Msort1(int l,int r){
	int mid=(l+r)>>1ll;
	int j=l,k=mid+1;
	for (register int i=l;i<=r;++i){
		if (p[i].id<=mid) t[j++]=p[i];
		else t[k++]=p[i];
	}
	for (register int i=l;i<=r;++i) p[i]=t[i];
}
inline void Msort2(int l,int r){
	int mid=(l+r)>>1ll;
	int j=l,k=mid+1;
	for (register int i=l;i<=r;++i){
		if (j<=mid&&(k>r||cmp2(p[j],p[k]))) t[i]=p[j++];
		else t[i]=p[k++];
	}
	for (register int i=l;i<=r;++i) p[i]=t[i];
}
void CDQ(int l,int r){
	if (l==r){
		if (l!=1) ans[l]=max(ans[l],ans[l-1]);
		p[l].x=ans[l]/(A[l]*R[l]+B[l]);
		p[l].y=p[l].x*R[l];
		return ;
	}
	int mid=(l+r)/2;
	Msort1(l,r);
	CDQ(l,mid);
	int head=1,rear=0;
	for (register int i=l;i<=mid;++i){
		while (head<rear&&Slope(q[rear]-q[rear-1])<=Slope(p[i]-q[rear-1])) rear--;
		q[++rear]=p[i];
	}
	for (register int k=mid+1;k<=r;++k){
		int i=p[k].id;
		while (head<rear&&Slope(q[head+1]-q[head])>=(-B[i]/A[i])) head++;
		int j=q[head].id;
		ans[i]=max(ans[i],(ans[j]/(A[j]*R[j]+B[j]))*(A[i]*R[j]+B[i]));
	}
	CDQ(mid+1,r);
	Msort2(l,r);
}
#undef int
int main(){
#define int long long
	N=read(),S=(double)read();
	for (register int i=1;i<=N;++i){
		cin>>A[i]>>B[i]>>R[i];
		p[i].id=i;
	}
	ans[1]=S;
	sort(p+1,p+1+N,cmp1);
	CDQ(1,N);
	printf("%.3f\n",ans[N]);
}

例题4：

Machine Works

题意：

一个公司获得了一个厂房n(10^5)天的使用权
和一笔启动资金C(10^9)，准备在n天里租借机器生产来获得收益
可以租借的机器有M(10^5)个，每个机器有四个值，D,P,R,G (D<=n, P,R,G都是10^9)
表明你可以再第D天花费P费用（首先手里必须有那么多钱）
租借这个机器，从D+1天开始该机器每天产生G的收益，在你不需要机器时
可以卖掉这个机器，一次获得R的钱

显然可以列出dp方程：$f[i] = max(f[j] - P[i] + R[j] + G[j] * (D[i] - D[j] - 1)) $

注意到f[j]>=0f[j]>=0才能转移。

代表你在第jj天买进机器，然后后面D[i]−D[j]−1D[i]-D[j]-1天使用机器，最后第ii 天卖出。

斜率优化比较好推，我就不推了。

代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
#include <bits/stdc++.h>
#define MAXN 100005
#define int long long
#define INF 1e15
#define eps 1e-10
using namespace std;
inline int read(){
	int x=0,f=1;
	char ch=getchar();
	while (ch<'0'||ch>'9'){
		if (ch=='-') f=-1;
		ch=getchar();
	}
	while (ch>='0'&&ch<='9'){
		x=(x<<3)+(x<<1)+(ch^'0');
		ch=getchar();
	}
	return x*f;
}
struct Machine{
	int d,p,r,g;
}m[MAXN];
inline bool cmp(const Machine &A,const Machine &B){
	return A.d<B.d;
}
int n,c,d;
int ans[MAXN];
struct Point{
	double x,y;
	int id;
}p[MAXN],t[MAXN],q[MAXN];
inline bool cmp1(const Point &A,const Point &B){
	return m[A.id].d<m[B.id].d;
}
inline bool cmp2(const Point &A,const Point &B){
	if (A.x!=B.x) return A.x<B.x;
	else return A.y<B.y;
}
inline double operator * (const Point &A,const Point &B){
	return A.x*B.y-A.y*B.x;
}
inline Point operator - (const Point &A,const Point &B){
	return Point{A.x-B.x,A.y-B.y};
}
inline Point operator + (const Point &A,const Point &B){
	return Point{A.x+B.x,A.y+B.y};
}
inline double Slope(const Point &A,const Point &B){
	return (double)(A.y-B.y)/(double)(A.x-B.x);
}
inline void Msort1(int l,int r){
	int mid=(l+r)>>1;
	int j=l,k=mid+1;
	for (register int i=l;i<=r;++i){
		if (p[i].id<=mid) t[j++]=p[i];
		else t[k++]=p[i];
	}
	for (register int i=l;i<=r;++i) p[i]=t[i];
}
inline void Msort2(int l,int r){
	int mid=(l+r)>>1;
	int j=l,k=mid+1;
	for (register int i=l;i<=r;++i){
		if (j<=mid&&(k>r||cmp2(p[j],p[k]))) t[i]=p[j++];
		else t[i]=p[k++];
	}
	for (register int i=l;i<=r;++i) p[i]=t[i];
}
void CDQ(int l,int r){
	if (l==r){
		p[l].x=(double)m[l].g;
		p[l].y=(double)ans[l]+m[l].r-m[l].g*m[l].d-m[l].g;
		return ;
	}
	int mid=(l+r)>>1;
	Msort1(l,r);
	CDQ(l,mid);
	int head=1,rear=0;
	for (register int i=l;i<=mid;++i){
		if (ans[p[i].id]<0) continue;
		while (head<rear&&(q[rear]-q[rear-1])*(p[i]-q[rear-1])>=0) rear--;
		q[++rear]=p[i];
	}
	for (register int k=mid+1;k<=r;++k){
		int i=p[k].id;
		while (head<rear&&(q[head+1]-q[head])*(Point{1,-m[i].d,1926})<=0) head++;
		int j=q[head].id;
		ans[i]=max(ans[i],ans[j]+m[j].r+m[j].g*(m[i].d-m[j].d-1)-m[i].p);
	}
	CDQ(mid+1,r);
	Msort2(l,r);
}
#undef int
int main(){
#define int long long
	int Case=0;
	while (scanf("%lld%lld%lld",&n,&c,&d)!=EOF&&(n!=0&&c!=0&&d!=0)){
		for (register int i=1;i<=n;++i){
			m[i].d=read(),m[i].p=read(),m[i].r=read(),m[i].g=read();
			p[i].id=i;
		}
		sort(m+1,m+1+n,cmp);
		sort(p+1,p+1+n,cmp1);
		for (register int i=1;i<=n;++i) ans[i]=-INF;
		ans[0]=c;
		CDQ(0,n);
		long long ret=c;
		for (register int i=1;i<=n;++i){
			if (ans[i]>=0) ret=max(ret,(long long)(ans[i]+m[i].r+(d-m[i].d)*m[i].g));
		}
		printf("Case %d: %lld\n",++Case,(long long)ret);
	}
}
相关文章
斜率优化学习笔记基本概念

博客内容遵循 署名-非商业性使用-相同方式共享 4.0 国际 (CC BY-NC-SA 4.0) 协议

本文永久链接是：https://gaisaiyuno.github.io/archives/b7a912cc.html

更新于：2022年6月14日

CDQ分治

斜率优化

评论