前言
这种方法比传统斜率优化更快,更准,更狠。
凸包优化
一切形如
d
p
[
i
]
=
min
/
max
{
f
1
(
j
)
⋅
g
1
(
i
)
+
f
2
(
j
)
}
+
g
2
(
i
)
dp[i]=\min/\max\{f_1(j) \cdot g_1(i) + f_2(j)\} + g_2(i)
dp[i]=min/max{f1(j)⋅g1(i)+f2(j)}+g2(i)的转移方程,都可以凸包优化。
其中,
f
f
f为关于
j
j
j的函数,
g
g
g为关于
i
i
i的函数。
例如
d
p
[
i
]
=
min
{
−
2
h
j
⋅
h
i
+
h
j
2
+
d
p
[
j
]
}
+
a
i
+
h
i
2
dp[i] = \min\{-2h_j \cdot h_i + {h_j}^2 + dp[j]\} + a_i + {h_i}^2
dp[i]=min{−2hj⋅hi+hj2+dp[j]}+ai+hi2(这里面,
f
1
(
j
)
=
−
2
h
j
f_1(j) = -2h_j
f1(j)=−2hj,
f
2
(
j
)
=
h
j
2
+
d
p
[
j
]
f_2(j) = {h_j}^2 + dp[j]
f2(j)=hj2+dp[j],
g
1
(
i
)
=
h
i
g_1(i) = h_i
g1(i)=hi,
g
2
(
i
)
=
a
i
+
h
i
2
g_2(i) = a_i + {h_i}^2
g2(i)=ai+hi2)
我们接下来口胡 d p [ i ] = max { f 1 ( j ) ⋅ g 1 ( i ) + f 2 ( j ) } + g 2 ( i ) dp[i]=\max\{f_1(j) \cdot g_1(i) + f_2(j)\} + g_2(i) dp[i]=max{f1(j)⋅g1(i)+f2(j)}+g2(i)的情况。
很简单。
第一步
定义一个关于 g 1 ( i ) g_1(i) g1(i)和 j j j的二元函数: L j ( g 1 ( i ) ) = f 1 ( j ) ⋅ g 1 ( i ) + f 2 ( j ) L_j\left(g_1(i)\right)=f_1(j) \cdot g_1(i) + f_2(j) Lj(g1(i))=f1(j)⋅g1(i)+f2(j)为什么叫 L j L_j Lj呢,因为这是一条直线,这条直线的斜率为 f 1 ( j ) f_1(j) f1(j),纵截距为 f 2 ( j ) f_2(j) f2(j)。
第二步
d p [ i ] = max { L j ( g 1 ( i ) ) } + g 2 ( i ) dp[i]=\max\{L_j(g_1(i))\} + g_2(i) dp[i]=max{Lj(g1(i))}+g2(i)也就是说,我们只需要找直线 x = g 1 ( i ) x = g_1(i) x=g1(i)与所有 L j L_j Lj的交点中纵坐标最大的那个。
最后一步
用个李超线段树即可。
但是,在大多数题你都会发现,
f
1
f_1
f1和
g
1
g_1
g1有单调性。
否则,用李超线段树或CDQ或平衡树什么的即可。
那么我接下来讲
f
1
f_1
f1单调减,
g
1
g_1
g1单调增的情况吧。
再说一遍,很简单。(你发现我们没有进行任何计算)
现在要计算
d
p
[
i
]
dp[i]
dp[i],则我们可以做到:此时已经按顺序把所有
L
j
(
1
≤
j
<
i
)
L_j(1\leq j <i)
Lj(1≤j<i)放进了一个双端队列
Q
Q
Q,呈这个样子(
1
1
1指
L
Q
[
H
e
a
d
]
L_{Q[Head]}
LQ[Head],
2
2
2指
L
Q
[
H
e
a
d
+
1
]
L_{Q[Head + 1]}
LQ[Head+1],以此类推):
加粗的地方是这个直线的“贡献”,但有些直线没有贡献,例如下图中的黑线:
基于归纳的思想,我们可以假设此时队列中没有这种线
(
∗
)
(*)
(∗),然后在该次DP后维护这样一个双端队列
Q
Q
Q。
一个显然的结论是:由于
g
1
(
i
)
g_1(i)
g1(i)单增,那么如果
g
1
(
i
)
g_1(i)
g1(i)到了这个地方,蓝线就没用了:
所以,不断比较
L
Q
[
H
e
a
d
]
(
g
1
(
i
)
)
L_{Q[Head]}(g_1(i))
LQ[Head](g1(i))和
L
Q
[
H
e
a
d
+
1
]
(
g
1
(
i
)
)
L_{Q[Head + 1]}(g1(i))
LQ[Head+1](g1(i)),来看
L
Q
[
H
e
a
d
]
L_{Q[Head]}
LQ[Head]有没有存在的必要,类似传统斜率优化。
然后,考虑加入当前直线
L
i
L_i
Li(下图中的黑色),如果是这样的,那么绿线就没有用了(
3
3
3指
L
Q
[
T
a
i
l
]
L_{Q[Tail]}
LQ[Tail],
2
2
2指
L
Q
[
T
a
i
l
−
1
]
L_{Q[Tail - 1]}
LQ[Tail−1],以此类推):
这个问题的刻画也很好想到:是比较
L
i
L_i
Li与
L
Q
[
T
a
i
l
]
L_{Q[Tail]}
LQ[Tail]的交点 和
L
Q
[
T
a
i
l
]
L_{Q[Tail]}
LQ[Tail]与
L
Q
[
T
a
i
l
−
1
]
L_{Q[Tail - 1]}
LQ[Tail−1]的交点 的横坐标。下图中,若
x
A
<
x
B
x_A<x_B
xA<xB,那
L
Q
[
T
a
i
l
]
L_{Q[Tail]}
LQ[Tail]就没用了:
于是这样就能做到
(
∗
)
(*)
(∗)了,也是类似于传统斜率优化。
说完了,看例题代码有惊♂喜。
例一
Kalila and Dimna in the Logging Industry
转移方程
不用看题,直接看转移方程即可: d p [ i ] = min 1 ≤ j < i { d p [ j ] + b j ⋅ a i } dp[i] = \min\limits_{1 \leq j < i}\{dp[j] + b_j \cdot a_i\} dp[i]=1≤j<imin{dp[j]+bj⋅ai}其中 a i a_i ai递增, b i b_i bi递减。
凸包优化
f 1 ( j ) = b j f_1(j)=b_j f1(j)=bj, g 1 ( i ) = a i g_1(i)=a_i g1(i)=ai, f 2 ( j ) = d p [ j ] f_2(j)=dp[j] f2(j)=dp[j], g 2 ( i ) = 0 g_2(i)=0 g2(i)=0,其中 f 1 f_1 f1单减, g 1 g_1 g1单增,跟上面讲的情况一模一样。
代码
#include <algorithm>
#include <cstdio>
#include <cstring>
typedef long long LL;
const int MAXN = 100000;
const LL INF = 1ll << 60;
int N;
LL A[MAXN + 5], B[MAXN + 5];
LL Dp[MAXN + 5];
struct Line {
LL k, b;
Line() { }
Line(LL _k, LL _b) { k = _k, b = _b; }
LL Calc(int x) { return k * x + b; } // 算函数值
double Ints(Line other) { // 求两直线交点的横坐标
return (double)(other.b - b) / (k - other.k);
}
}Q[MAXN + 5];
int Head, Tail;
int main() {
scanf("%d", &N);
for (int i = 1; i <= N; i++)
scanf("%lld", A + i);
for (int i = 1; i <= N; i++)
scanf("%lld", B + i);
Q[Head = Tail = 1] = Line(B[1], 0); // 边界注意一下即可
for (int i = 2; i <= N; i++) {
int x = A[i];
while (Tail - Head + 1 >= 2 && Q[Head].Calc(x) >= Q[Head + 1].Calc(x))
Head++;
Dp[i] = Q[Head].Calc(x); // 找到x=A[i]处的最低点
Line cur(B[i], Dp[i]);
while (Tail - Head + 1 >= 2 && Q[Tail].Ints(cur) <= Q[Tail].Ints(Q[Tail - 1]))
Tail--;
Q[++Tail] = cur; // 加入Li
}
printf("%lld\n", Dp[N]);
return 0;
}
例二
题目大意
你想打开
z
z
z个椰子吃,你的沙比队友给你准备了
n
n
n个椰子,每个椰子的坚硬♂程度不同,第
i
i
i个椰子的坚硬♂程度是
a
i
a_i
ai,表示它要被敲
a
i
a_i
ai下才能被打开(不一定要连续敲)。 你不知道椰子的顺序。 请问至少要敲多少下才能打开最少
z
z
z个椰子。
有必要看一下样例:
Input
2
2 1
50 55
2 1
40 100
Output
55
80
第一个:抓一个直接敲55下,不管怎么样都能敲开;
第二个:抓一个,先敲40下,如果没开,就拿另一个敲40下,至少能得到1个椰子。
转移方程
我都没看出来是个DP。
先排个序,然后先考虑怎么敲开一个椰子:
记阴影矩形的面积为
S
i
S_i
Si,如果我们想撬开1个椰子,那敲
min
{
S
i
}
\min\{S_i\}
min{Si}下就行了,因为对于任意一种
a
i
×
(
n
−
i
+
1
)
a_i\times(n-i+1)
ai×(n−i+1)下的方案,必定能敲出一个椰子:先随便找个椰子敲
a
i
a_i
ai下,如果没打开,就换一个没敲过的再敲,重复此操作,脸再黑也就是把阴影部分倒着敲完,那也能把第
i
i
i个敲开。
接下来考虑,如果我们想敲开两个椰子,答案是
min
i
<
j
{
S
i
∪
S
j
}
\min\limits_{i<j}\{S_i\cup S_j\}
i<jmin{Si∪Sj}。
考虑你是一个黑人的情况:先敲了
S
i
S_i
Si下才敲开一个椰子,那你的椰子变成了这样:
然后,你肯定知道哪些是敲过的,你就在敲过的那些里面敲
S
j
S_j
Sj下,就又打开了一个椰子。
于是问题转变为在矩形里面找面积最小的,含 z z z级的阶梯的阶梯形(我是倒着来的): d p [ i ] [ j ] = min k > j { d p [ i − 1 ] [ k ] + a j ⋅ ( k − j ) } dp[i][j] = \min\limits_{k>j}\{dp[i - 1][k] + a_j\cdot(k-j)\} dp[i][j]=k>jmin{dp[i−1][k]+aj⋅(k−j)}
凸包优化
d
p
[
i
]
[
j
]
=
min
k
>
j
{
d
p
[
i
−
1
]
[
k
]
+
k
⋅
a
j
}
−
a
j
⋅
j
dp[i][j] = \min\limits_{k>j}\{dp[i - 1][k] + k\cdot a_j\}-a_j\cdot j
dp[i][j]=k>jmin{dp[i−1][k]+k⋅aj}−aj⋅j
f
1
(
k
)
=
k
f_1(k)=k
f1(k)=k,
f
2
(
k
)
=
d
p
[
i
−
1
]
[
k
]
f_2(k)=dp[i - 1][k]
f2(k)=dp[i−1][k],
g
1
(
j
)
=
a
j
g_1(j) = a_j
g1(j)=aj,
g
2
(
j
)
=
a
j
⋅
j
g_2(j) = a_j\cdot j
g2(j)=aj⋅j,注意
i
i
i跟凸包优化无关,是
j
,
k
j,k
j,k参与凸包优化。
由于我倒着来的,所以
f
1
f_1
f1单减,
g
1
g_1
g1单减,然后就简单了。
代码
#include <algorithm>
#include <cstdio>
#include <cstring>
typedef long long LL;
const int MAXN = 1000;
int N, Z; LL H[MAXN + 5];
LL Dp[MAXN + 5][MAXN + 5];
struct Line {
LL k, b;
Line() { }
Line(LL x, LL y) { k = x, b = y; }
LL Calc(int x) {
return k * x + b;
}
double Ints(Line other) {
return (double)(b - other.b) / (other.k - k);
}
}Q[MAXN + 5];
int Head, Tail;
/*
1
3 2
1 8 10
*/
int main() {
int T; scanf("%d", &T);
while (T--) {
scanf("%d%d", &N, &Z);
for (int i = 1; i <= N; i++)
scanf("%lld", &H[i]);
std::sort(H + 1, H + 1 + N);
for (int i = 1; i <= N; i++)
Dp[1][i] = (N - i + 1) * H[i];
for (int i = 2; i <= Z; i++) {
Q[Head = Tail = 1] = Line(N - i + 2, Dp[i - 1][N - i + 2]);
for (int j = N - i + 1; j >= 1; j--) { // 注意边界
int x = H[j];
while (Tail - Head + 1 >= 2 && Q[Tail].Calc(x) >= Q[Tail - 1].Calc(x))
Tail--;
Dp[i][j] = Q[Tail].Calc(x) - H[j] * j;
Line cur(j, Dp[i - 1][j]); // 当前层是加上一层的直线 通过转移方程就能看出来
while (Tail - Head + 1 >= 2 && Q[Tail].Ints(cur) <= Q[Tail].Ints(Q[Tail - 1]))
Tail--;
Q[++Tail] = cur;
}
}
LL Ans = 1ll << 60;
for (int i = 1; i <= N - Z + 1; i++)
Ans = std::min(Ans, Dp[Z][i]);
printf("%lld\n", Ans);
}
return 0;
}
李超线段树
如果
f
1
f_1
f1和
g
1
g_1
g1没有单调性,我们就不能用双端队列维护了。
李超线段树的作用很简单:维护一些一次函数(直线 / 线段),支持插入和查询,查询时可以找到当前横坐标下最大 / 最小的函数值。
完美解决几乎所有凸包优化。
代码只有40行。
思想
它每个区间记录的是该区间中点处的最大函数值对应的函数 M a x i Max_i Maxi。
插入
插入直线 c u r cur cur的过程如下:
- c u r cur cur在这个区间上完全覆盖了 M a x i Max_i Maxi:将 M a x i Max_i Maxi变成 c u r cur cur,返回(没有懒标记,不用再改儿子,看查询的过程就知道了);
- 如果该区间中点处 M a x i ( m i d ) < c u r ( m i d ) Max_i(mid)<cur(mid) Maxi(mid)<cur(mid),则交换 M a x i Max_i Maxi和 c u r cur cur,保证 M a x i Max_i Maxi的意义正确;
- 现在的 c u r cur cur会对交点所在子树产生贡献(下图中,右子树的橙色段需要修改),因此递归下去:
查询
比较简单,递归得到下层的答案,跟自己这层比(因此不用插入和查询都可以不用懒标记)即可。
代码
见例题,有惊♂喜。
例三
[JSOI2008]Blue Mary开公司
这是一道版题。
代码
#include <algorithm>
#include <cstdio>
#include <cstring>
const int MAXT = 100000;
const int MAXX = 50000;
const double INF = 1e9;
struct LiChao_Tree {
#define lch (i << 1)
#define rch (i << 1 | 1)
struct Line {
double k, b;
inline double Calc(int x) {
return k * x + b;
}
}Max[MAXT + 5];
inline bool Cover(Line Low, Line High, int x) { // 判断x处Hight否覆盖了Low
return Low.Calc(x - 1) <= High.Calc(x - 1);
}
void Insert(int i, int l, int r, Line cur) {
if (Cover(Max[i], cur, l) && Cover(Max[i], cur, r)) {
Max[i] = cur;
return;
}
if (l == r)
return;
int mid = (l + r) >> 1;
if (Cover(Max[i], cur, mid))
std::swap(Max[i], cur);
if (Cover(Max[i], cur, l))
Insert(lch, l, mid, cur);
if (Cover(Max[i], cur, r))
Insert(rch, mid + 1, r, cur);
}
double Query(int i, int l, int r, int x) {
double tmp = -INF;
int mid = (l + r) >> 1;
if (x < mid)
tmp = Query(lch, l, mid, x);
if (x > mid)
tmp = Query(rch, mid + 1, r, x);
return std::max(tmp, Max[i].Calc(x - 1));
}
}Tree;
int main() {
int T, X; scanf("%d", &T);
while (T--) {
char opt[20];
scanf("%s", opt);
if (opt[0] == 'P') {
LiChao_Tree::Line tmp;
scanf("%lf%lf", &tmp.b, &tmp.k);
Tree.Insert(1, 1, MAXX, tmp);
}
else {
scanf("%d", &X);
printf("%d\n", int(Tree.Query(1, 1, MAXX, X) / 100));
}
}
return 0;
}
例四
转移方程
d p [ i ] = min j < i 且 p j < p i { d p [ j ] + ( h i − h j ) 2 } + a i dp[i]=\min\limits_{_{j<i\text{且}p_j<p_i}}\{dp[j]+(h_i-h_j)^2\}+a_i dp[i]=j<i且pj<pimin{dp[j]+(hi−hj)2}+ai其中 p p p不单调, h h h不单调, a a a不单调。
怎么做
看到这个题,什么都不单调,还尼玛有转移限制???
不可做,溜了。
正解:树状数组套李超树维护凸包。
树状数组中,每个结点是一个李超树,维护对应区间的凸包。查询的时候,从
p
i
p_i
pi用lowbit
减到
0
0
0,根据树状数组的性质,访问到的恰好就是
d
p
[
i
]
dp[i]
dp[i]的所有转移直线,统计最大的函数值即可。(其实树状数组很大的一个用处就是处理偏序问题,一定程度上可以替代CDQ分治)
代码
#include <algorithm>
#include <cstdio>
#include <cstring>
typedef long long LL;
const int MAXN = 300000;
const int MAXL = 600000;
const LL INF = 1ll << 60;
struct Line {
LL k, b;
Line() { k = 0, b = INF; }
Line(LL _k, LL _b) { k = _k, b = _b; }
LL Calc(int x) { return k * x + b; }
double Ints(Line other) {
return (double)(other.b - b) / (k - other.k);
}
};
struct LiChao_Tree {
#define lch (Child[i][0])
#define rch (Child[i][1])
Line Min[MAXN * 20 + 5];
int NodeCnt;
int Child[MAXN * 20 + 5][2];
inline bool Cover(Line Low, Line High, int x) {
return Low.Calc(x) <= High.Calc(x);
}
void Insert(int &i, int l, int r, Line cur) {
if (!i)
i = ++NodeCnt;
if (Cover(cur, Min[i], l) && Cover(cur, Min[i], r)) {
Min[i] = cur;
return;
}
if (l == r)
return;
int mid = (l + r) >> 1;
if (Cover(cur, Min[i], mid))
std::swap(Min[i], cur);
if (Cover(cur, Min[i], l))
Insert(lch, l, mid, cur);
if (Cover(cur, Min[i], r))
Insert(rch, mid + 1, r, cur);
}
LL Query(int i, int l, int r, int x) {
LL tmp = INF;
int mid = (l + r) >> 1;
if (x < mid)
tmp = Query(lch, l, mid, x);
if (x > mid)
tmp = Query(rch, mid + 1, r, x);
return std::min(tmp, Min[i].Calc(x));
}
#undef lch
#undef rch
}Tree;
struct BIT {
#define lowbit(x) ((x) & (-(x)))
int Root[MAXN + 5];
void Update(int p, Line l) {
for (int i = p; i <= MAXN; i += lowbit(i))
Tree.Insert(Root[i], 1, MAXL, l);
}
LL GetMin(int p, int x) {
LL ret = INF;
for (int i = p; i > 0 ; i -= lowbit(i))
ret = std::min(ret, Tree.Query(Root[i], 1, MAXL, x));
return ret;
}
#undef lowbit
}CHT;
int N, P[MAXN + 5];
LL A[MAXN + 5], H[MAXN + 5];
LL Dp[MAXN + 5];
int main() {
scanf("%d", &N);
for (int i = 1; i <= N; i++)
scanf("%d", &P[i]);
for (int i = 1; i <= N; i++)
scanf("%lld", &A[i]);
for (int i = 1; i <= N; i++)
scanf("%lld", &H[i]);
CHT.Update(P[1], Line(-2 * H[1], A[1] + H[1] * H[1]));
for (int i = 2; i <= N; i++) {
Dp[i] = CHT.GetMin(P[i], H[i]) + A[i] + H[i] * H[i];
CHT.Update(P[i], Line(-2 * H[i], Dp[i] + H[i] * H[i]));
}
printf("%lld", Dp[N]);
return 0;
}