参考资料:
zzq yyds,算法导论,董克凡、邹逍遥的论文
由于博主实力过菜,所以大部分都是参考的,其中算法导论的讲解最为详尽.
但本文依然希望能用尽量精简的说明来让读者理解.
定义
已知一组实数
a
1
,
a
2
.
.
.
.
a
n
a_1,a_2....a_n
a1,a2....an,一组变量
x
1
,
x
2
.
.
.
x
n
x_1,x_2...x_n
x1,x2...xn,
在这些变量上的 线性函数 定义为
f
(
x
1
,
x
2
,
.
.
.
,
x
n
)
=
∑
i
=
1
n
a
i
x
i
f(x_1,x_2,...,x_n)=\sum_{i=1}^n a_ix_i
f(x1,x2,...,xn)=∑i=1naixi.
如果用向量 A , X A,X A,X表示两个数列.则可简记为 f ( X ) = A T X f(X)=A^T X f(X)=ATX.
定义等式 f ( X ) = b f(X)=b f(X)=b,不等式 f ( X ) ≤ b o r f ( X ) ≥ b f(X)\le b ~or~ f(X)\ge b f(X)≤b or f(X)≥b为线性约束.
则线性规划定义为在满足有限线性约束的情况下最优化一个线性函数(目标函数).
线性规划有两种形式: 标准型,松弛型.
- 标准型的优点是简洁明了,缺点在于其在移项的时候可能需要符号的变化,对于程序实现很不方便.
松弛型给每一个线性约束都定义了一个变量 x i + n x_{i+n} xi+n,容易发现 x i + n ≥ 0 x_{i+n}\ge 0 xi+n≥0和原来的约束等价.
标准型中的所有约束均为符号是 ≤ \le ≤的线性不等式,松弛型的所有约束则表示为线性等式. - 在松弛型中,我们定义
x
i
+
n
,
i
=
1
,
2...
m
x_{i+n},i=1,2...m
xi+n,i=1,2...m为基变量,
x
i
,
i
=
1...
n
x_i,i=1...n
xi,i=1...n为非基变量.
可以发现松弛型的系数隐含着一组解:
令所有非基变量的值取 0,基变量 x i + n = b i x_{i+n}=b_i xi+n=bi,我们称这样的解为基本解.
∀ b i ≥ 0 \forall b_i\ge 0 ∀bi≥0,则容易看出这组解是可行的,称为基本可行解.
对于后面的算法部分我们只研究基本解.
线性规划的性质
- 任意一种线性规划都可以转化为标准型的形式:
对于 f ( X ) ≤ b f(X)\le b f(X)≤b,直接加入.
对于 f ( X ) ≥ b f(X)\ge b f(X)≥b,加入 − f ( X ) ≤ − b -f(X)\le -b −f(X)≤−b.
对于 f ( X ) = b f(X)=b f(X)=b,可以拆成 f ( X ) ≤ b a n d f ( X ) ≥ b f(X)\le b ~and ~f(X)\ge b f(X)≤b and f(X)≥b
而对于最小化 f ( X ) f(X) f(X)的问题,我们也可以转化为最大化 − f ( X ) -f(X) −f(X). - 从几何角度观察:
如果我们把一个解看作一个 n n n维的点 ( x 1 , x 2 . . . x n ) (x_1,x_2...x_n) (x1,x2...xn),
那么对于每一个线性约束,我们可以发现它对应一个超平面(自由度为 n − 1 n-1 n−1),
而解空间就可以表示为超平面的交,这个交是一个凸形区域(感性理解为区域内任意两点间的连线上的点在区域内),而一个凸形区域内 f ( X ) f(X) f(X)的局部最优解是全局最优解,所以我们只要从凸形区域的一个顶点走到另一个顶点即可 ↓ \downarrow ↓.
(论文说用归纳证明,但我不是很懂.如果有比较初等的证明请指教)
一个简单的例子:
对于上面这个有两个变量的线性规划,我们要最大化 x 1 + x 2 x_1+x_2 x1+x2,由于约束的交形成了一个凸包,所以顶点一定是最优决策.
单纯形法
单纯形法是求解线性规划的有利工具.
虽然线性规划有多项式级别的算法,但是指数过高,在 O I OI OI中不实用.
单纯形法需要我们把线性规划转为松弛型.
之后,算法有一个重要的操作
p
i
v
o
t
(
l
,
e
)
pivot(l,e)
pivot(l,e)(转轴),即替换
x
l
+
n
,
x
e
x_{l+n},x_e
xl+n,xe使得
x
l
+
n
x_{l+n}
xl+n变为非基变量,
x
e
x_e
xe变为基变量.
设基变量的下标集合为
B
B
B,非基变量的下标集合为
N
N
N.
x
l
+
n
=
b
l
−
∑
j
∈
N
a
l
,
j
x
j
x_{l+n}=b_l-\sum_{j\in N} a_{l,j}x_j
xl+n=bl−j∈N∑al,jxj
x l + n = b l − ∑ j ∈ N , j ≠ e a l , j x j − a l , e x e x_{l+n}=b_l-\sum_{j\in N,j\ne e} a_{l,j}x_j-a_{l,e}x_e xl+n=bl−j∈N,j=e∑al,jxj−al,exe
a l , e x e = b l − ∑ j ∈ N , j ≠ e a l , j x j − x l + n a_{l,e}x_e=b_l-\sum_{j\in N,j\ne e} a_{l,j}x_j-x_{l+n} al,exe=bl−j∈N,j=e∑al,jxj−xl+n
a l , e x e = b l − ∑ j ∈ N , j ≠ e a l , j x j − x l + n a_{l,e}x_e=b_l-\sum_{j\in N,j\ne e} a_{l,j}x_j-x_{l+n} al,exe=bl−j∈N,j=e∑al,jxj−xl+n
x e = b l − ∑ j ∈ N , j ≠ e a l , j x j − x l + n a l , e x_e=\dfrac{b_l-\sum_{j\in N,j\ne e} a_{l,j}x_j-x_{l+n}}{a_{l,e}} xe=al,ebl−∑j∈N,j=eal,jxj−xl+n
然后我们把所有约束和目标函数中的 x e x_e xe换成上式,复杂度为 O ( n m ) O(nm) O(nm).
假设存在一组解
(
0
,
0
,
0..
,
0
,
b
1
,
b
2
.
.
.
b
m
)
(0,0,0..,0,b_1,b_2...b_m)
(0,0,0..,0,b1,b2...bm),即
∀
b
i
≥
0
\forall b_i\ge 0
∀bi≥0(忽略
i
n
i
t
i
a
l
i
z
a
t
i
o
n
initialization
initialization),那么我们执行以下操作:
若
∀
i
∈
B
a
i
,
e
≤
0
\forall_{i\in B} a_{i,e}\le 0
∀i∈Bai,e≤0,那么我们令
x
i
=
∞
x_i=\infty
xi=∞显然和所有约束不冲突且目标函数取得
∞
\infty
∞,故返回
"
u
n
b
o
u
n
d
e
d
"
"unbounded"
"unbounded"(无界).
否则,执行转轴后依然得到一组基本解(
b
i
′
=
b
i
−
a
i
,
e
x
e
≥
0
)
b_i'=b_i-a_{i,e}x_e\ge 0)
bi′=bi−ai,exe≥0),
同时由于我们
x
e
x_e
xe增大到对于所有约束最紧的值,所以目标函数不会变劣.
但是存在一种情况即
b
l
=
0
b_l=0
bl=0使得目标函数不变.
这个时候我们应有
B
l
a
n
d
Bland
Bland规则:对于同值取最小的下标,就可以使得这个确定性算法不会循环(证明不会,这里是引用算法导论的(算法导论也略去了证明)).(终止性)
几何意义:这里本质上是相当于沿着超平面交的边走到另一个顶点,然后局部最优的选择策略可以导向全局最优.
最优性代数证明:算法导论给出了一种构造对偶线性规划的方法,并构造对偶线性规划的一组合法解满足两个目标函数相等,从而得出结束循环时得到目标函数的最优值.
线性规划对偶定理:线性规划的目标函数最大值等于其对偶线性规划的目标函数最小值.
由于篇幅过长,感兴趣的读者自己去看证明吧~~
上述解法的大前提是
∀
b
i
≥
0
\forall b_i\ge 0
∀bi≥0,但是如果不满足咋办?
此时我们需要定义一个辅助线性规划来找到一组基本可行解.
容易发现存在一组可行解当且仅当辅助线性规划有一组解满足
x
0
=
0
x_0=0
x0=0.
我们要使得所有的
b
i
≥
0
b_i\ge 0
bi≥0,我们只要找到
b
i
b_i
bi最小的情况,然后
p
i
v
o
t
(
i
,
0
)
pivot(i,0)
pivot(i,0)即可.
然后我们执行伪代码中
3
∼
9
3\sim 9
3∼9的段即可.
如果
x
0
>
0
x_0>0
x0>0,则判断无解.
否则,若
x
0
x_0
x0为基变量,随便找到一个非基变量进行转轴.
最后移除
x
0
x_0
x0即可.
然后我们要把目标函数还原,原来的变量变为基变量的话我们把基变量展开,否则直接加系数即可.
至此,算法的大致内容讲完了.
最后再举一个小例子,来更好地理解线性规划的单纯形法:
辅助程序
{
最
小
化
x
1
+
x
2
+
x
3
+
x
4
满
足
约
束
:
−
2
x
1
+
8
x
2
+
0
x
3
+
10
x
4
≥
50
5
x
1
+
2
x
2
+
0
x
3
+
0
x
4
≥
100
3
x
1
−
5
x
2
+
10
x
3
−
2
x
4
≥
25
x
1
,
x
2
,
x
3
,
x
4
≥
0
\begin{cases} 最小化 & x_1+x_2+x_3+x_4 \\ 满足约束:&-2x_1+8x_2+0x_3+10x_4\ge 50\\ &5x_1+2x_2+0x_3+0x_4\ge 100\\ &3x_1-5x_2+10x_3-2x_4\ge 25\\ &x_1,x_2,x_3,x_4\ge 0 \end{cases}
⎩⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎧最小化满足约束:x1+x2+x3+x4−2x1+8x2+0x3+10x4≥505x1+2x2+0x3+0x4≥1003x1−5x2+10x3−2x4≥25x1,x2,x3,x4≥0
先转化为松弛型:
{
最
大
化
−
x
1
−
x
2
−
x
3
−
x
4
x
5
=
−
50
−
2
x
1
+
8
x
2
+
0
x
3
+
10
x
4
x
6
=
−
100
+
5
x
1
+
2
x
2
+
0
x
3
+
0
x
4
x
7
=
−
25
+
3
x
1
−
5
x
2
+
10
x
3
−
2
x
4
x
1
,
x
2
,
x
3
,
x
4
.
.
.
x
7
≥
0
\begin{cases} 最大化 & -x_1-x_2-x_3-x_4 \\ &x_5=-50-2x_1+8x_2+0x_3+10x_4\\ &x_6=- 100+5x_1+2x_2+0x_3+0x_4\\ &x_7=- 25+3x_1-5x_2+10x_3-2x_4\\ &x_1,x_2,x_3,x_4...x_7\ge 0 \end{cases}
⎩⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎧最大化−x1−x2−x3−x4x5=−50−2x1+8x2+0x3+10x4x6=−100+5x1+2x2+0x3+0x4x7=−25+3x1−5x2+10x3−2x4x1,x2,x3,x4...x7≥0
构造辅助线性规划:
{
最
大
化
−
x
0
x
5
=
−
50
−
2
x
1
+
8
x
2
+
0
x
3
+
10
x
4
+
x
0
x
6
=
−
100
+
5
x
1
+
2
x
2
+
0
x
3
+
0
x
4
+
x
0
x
7
=
−
25
+
3
x
1
−
5
x
2
+
10
x
3
−
2
x
4
+
x
0
x
1
,
x
2
,
x
3
,
x
4
.
.
.
x
7
≥
0
\begin{cases} 最大化 &-x_0 \\ &x_5=-50-2x_1+8x_2+0x_3+10x_4+x_0\\ &x_6=- 100+5x_1+2x_2+0x_3+0x_4+x_0\\ &x_7=- 25+3x_1-5x_2+10x_3-2x_4+x_0\\ &x_1,x_2,x_3,x_4...x_7\ge 0 \end{cases}
⎩⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎧最大化−x0x5=−50−2x1+8x2+0x3+10x4+x0x6=−100+5x1+2x2+0x3+0x4+x0x7=−25+3x1−5x2+10x3−2x4+x0x1,x2,x3,x4...x7≥0
x
0
=
100
−
5
x
1
−
2
x
2
+
0
x
3
+
0
x
4
x0=100-5x_1-2x_2+0x_3+0x_4
x0=100−5x1−2x2+0x3+0x4
具体实现(97分是因为卡精度了…).
然而UOJ有一个初始化部分更加简单的实现worse.
以下为鄙人参考后写出的:
#include <bits/stdc++.h>
#define rep(i, a, b) for (int i = a; i <= b; i++)
#define REP(i, a, b) for (int i = b; i >= a; i--)
#define debug(x) cerr << #x << ' ' << '=' << ' ' << x << endl
using namespace std;
typedef double db;
const int N = 44;
const db eps=1e-8,inf=1e15;
int n,m,type,id[N],tp[N];
db a[N][N];
void pivot(int r,int c) {//r为基变量,c为非基变量.现在互换.为了使得目标函数更大.
swap(id[r+n],id[c]);
db t=-a[r][c]; a[r][c]=-1; rep(i,0,n) a[r][i] /= t;
rep(i,0,m) if(a[i][c]&&r^i) { t=a[i][c]; a[i][c]=0; rep(j,0,n) a[i][j] += t * a[r][j]; }
}
void solve() {
scanf("%d %d %d",&n,&m,&type);
rep(i,1,n) scanf("%lf",&a[0][i]);
rep(i,1,m) {rep(j,1,n) scanf("%lf",&a[i][j]),a[i][j] *= -1; scanf("%lf",&a[i][0]);}
db t;rep(i,1,n) id[i]=i;
while(1) {
int i=0,j=0; db w=-eps;
rep(k,1,m) if(a[k][0]<w) w=a[i=k][0];
if(!i) break;//初始化,使得所有b[i]>=0
rep(k,1,n) if(a[i][k]>eps) {j=k; break;}
if(!j) {puts("Infeasible"); return ;}//所有变量都是负系数则无解.
pivot(i,j);//转轴,使一个变量尽可能大.
}
while(1) {
int i=0,j=0; db w=eps;
rep(k,1,n) if(a[0][k]>w) w=a[0][j=k];//选择目标函数系数最大的变量来增大以增大函数.
if(!j) break; w=inf;
rep(k,1,m) if(a[k][j]<-eps && (t=-a[k][0]/a[k][j])<w) w=t,i=k;
//求一个能变成的最大值.
if(!i) {puts("Unbounded"); return ;}
pivot(i,j);
}
printf("%.9lf\n",a[0][0]);
rep(i,n+1,n+m) tp[id[i]]=i-n;
if(type) rep(i,1,n) printf("%.9lf ",tp[i]?a[tp[i]][0]:0);
}
int main() {
solve();
return 0;
}
空间复杂度为
O
(
n
m
)
O(nm)
O(nm),时间复杂度在数值为指数级别的时候会被卡成指数级,但是实际
O
I
OI
OI使用中期望
m
m
m次左右可以得出答案,所以可以像
d
i
n
i
c
dinic
dinic一样信仰跑.
重要性质:
当约束为线性等式时,
若系数矩阵为幺模矩阵,那么有最优整数解.
证明:运用克拉默法则,
x
i
=
det
A
i
det
A
x_i=\dfrac {\det A_i}{\det A}
xi=detAdetAi.
由于
det
A
=
±
1
\det A=\pm 1
detA=±1,所以基本解必为整数.
当约束为线性不等式时, A A A为全幺模矩阵那么极点为整数解.
例题
1430G - Yet Another DAG Problem
DAG的相关矩阵是全幺模矩阵,然后直接做就好.
我们把
x
i
x_i
xi当做
a
i
a_i
ai.
最小化
∑
i
=
1
m
w
i
∗
(
x
s
t
i
−
x
e
d
i
)
\sum_{i=1}^m w_i*(x_{st_i}-x_{ed_i})
∑i=1mwi∗(xsti−xedi).
满足约束:
x
i
≥
0
x
s
t
i
−
x
e
d
i
≥
1
x_i\ge 0\\x_{st_i}-x_{ed_i}\ge 1
xi≥0xsti−xedi≥1.
知道结论后就很裸了
#include<bits/stdc++.h>
#define rep(i,a,b) for(int i=a;i<=b;i++)
#define FOR(i,n) rep(i,1,n)
#define db double
using namespace std;
const int N=22,M=N*N;
const db eps=1e-8,inf=1e9;
int n,m,type=1,id[N+M],tp[N];
db a[M][N];
void pivot(int r,int c) {
swap(id[r+n],id[c]);
db t=-a[r][c]; a[r][c]=-1; rep(j,0,n) a[r][j] /= t;
rep(i,0,m) if(a[i][c] && i^r) {t=a[i][c]; a[i][c]=0; rep(j,0,n) a[i][j] += a[r][j]*t;}
}
void solve() {
db t;
FOR(i,n) id[i]=i;
while(1) {
int i=0,j=0; db w=-eps;
FOR(k,m) if(a[k][0]<w) w=a[i=k][0];
if(!i) break;
FOR(k,n) if(a[i][k]>eps) {j=k; break;}
if(!j) {puts("Infeasible"); return ;}
pivot(i,j);
}
while(1) {
int i=0,j=0; db w=eps;
FOR(k,n) if(a[0][k]>w) w=a[0][j=k];
if(!j) break; w=inf;
FOR(k,m) if(-a[k][j]>eps && (t=-a[k][0]/a[k][j]) < w) w=t,i=k;
if(!i) {puts("Unbounded"); return ;}
pivot(i,j);
}
a[0][0]=0;
rep(i,n+1,n+m) tp[id[i]]=i-n;
if(type) FOR(i,n) printf("%d ",(int)(a[tp[i]][0]+0.2));
puts("");
}
int main() {
scanf("%d %d",&n,&m);
FOR(i,m) {
int x,y,z; scanf("%d%d%d",&x,&y,&z);
a[0][x] -= z; a[0][y] += z;
a[i][x] = 1; a[i][y] = -1; a[i][0]=-1;
}
solve(); return 0;
}
P3980 [NOI2008]志愿者招募
设
x
j
x_j
xj为第
j
j
j类选择的人数.
则线性规划如下:
{
m
i
n
i
m
i
z
e
x
j
c
j
∑
A
i
,
j
x
j
≥
a
i
(
A
i
,
j
=
1
当
j
能
覆
盖
第
i
天
)
\begin{cases} minimize~~x_jc_j \\ \sum A_{i,j}x_j\ge a_i (A_{i,j}=1当j能覆盖第i天) \end{cases}
{minimize xjcj∑Ai,jxj≥ai(Ai,j=1当j能覆盖第i天)
这里我们可以直接转成松弛型也可以利用线性规划对偶定理.
直接转化:
#include<bits/stdc++.h>
#define rep(i,a,b) for(int i=a;i<=b;i++)
#define FOR(i,n) rep(i,1,n)
#define db double
using namespace std;
const int N=10010,M=1010;
const db eps=1e-8,inf=1e9;
int n,m,type,id[N+M],tp[N];
db a[M][N];
void pivot(int r,int c) {
swap(id[r+n],id[c]);
db t=-a[r][c]; a[r][c]=-1; rep(j,0,n) a[r][j] /= t;
rep(i,0,m) if(a[i][c] && i^r) {t=a[i][c]; a[i][c]=0; rep(j,0,n) a[i][j] += a[r][j]*t;}
}
void solve() {
db t;
FOR(i,n) id[i]=i;
while(1) {
int i=0,j=0; db w=-eps;
FOR(k,m) if(a[k][0]<w) w=a[i=k][0];
if(!i) break;
FOR(k,n) if(a[i][k]>eps) {j=k; break;}
if(!j) {puts("Infeasible"); return ;}
pivot(i,j);
}
while(1) {
int i=0,j=0; db w=eps;
FOR(k,n) if(a[0][k]>w) w=a[0][j=k];
if(!j) break; w=inf;
FOR(k,m) if(-a[k][j]>eps && (t=-a[k][0]/a[k][j]) < w) w=t,i=k;
if(!i) {puts("Unbounded"); return ;}
pivot(i,j);
}
printf("%d\n",(int)-a[0][0]);
}
void qr(int &x) {scanf("%d",&x);}
int main() {
qr(m); qr(n);
FOR(i,m) scanf("%lf",&a[i][0]),a[i][0] *= -1;
FOR(i,n) {
int x,y,z; qr(x); qr(y); qr(z); a[0][i]=-z;
rep(j,x,y) a[j][i]++;
}
solve(); return 0;
}
利用线性规划对偶定理,转化成:
{
m
a
x
i
m
i
z
e
y
j
a
j
∑
A
i
,
j
y
j
≤
c
i
(
A
j
,
i
=
1
当
j
能
覆
盖
第
i
天
)
\begin{cases} maximize~~y_ja_j \\ \sum A_{i,j}y_j\le c_i (A_{j,i}=1当j能覆盖第i天) \end{cases}
{maximize yjaj∑Ai,jyj≤ci(Aj,i=1当j能覆盖第i天)
#include<bits/stdc++.h>
#define rep(i,a,b) for(int i=a;i<=b;i++)
#define FOR(i,n) rep(i,1,n)
#define db double
using namespace std;
const int N=1010,M=10100;
const db eps=1e-8,inf=1e9;
int n,m,type,id[N+M],tp[N];
db a[M][N];
void pivot(int r,int c) {
swap(id[r+n],id[c]);
db t=-a[r][c]; a[r][c]=-1; rep(j,0,n) a[r][j] /= t;
rep(i,0,m) if(a[i][c] && i^r) {t=a[i][c]; a[i][c]=0; rep(j,0,n) a[i][j] += a[r][j]*t;}
}
void solve() {
db t;
FOR(i,n) id[i]=i;
while(1) {
int i=0,j=0; db w=eps;
FOR(k,n) if(a[0][k]>w) w=a[0][j=k];
if(!j) break; w=inf;
FOR(k,m) if(-a[k][j]>eps && (t=-a[k][0]/a[k][j]) < w) w=t,i=k;
if(!i) {puts("Unbounded"); return ;}
pivot(i,j);
}
printf("%d\n",(int)a[0][0]);
}
void qr(int &x) {scanf("%d",&x);}
int main() {
qr(n); qr(m);
FOR(i,n) scanf("%lf",&a[0][i]);
FOR(i,m) {
int x,y,z; qr(x); qr(y); qr(z);
a[i][0]=z; rep(j,x,y) a[i][j]--;
}
solve(); return 0;
}
可以发现,线性规划对偶定理在系数矩阵上的作用相当于进行一次转置然后系数取相反数,最后答案互为相反数.
当然这题也可以用费用流做.参考巨巨
我们可以发现这题和
k
k
k可重区间覆盖有点类似,但是现在是每个位置有一个下限.
我们容易想到连边
(
i
,
i
+
1
)
(i,i+1)
(i,i+1)来表示第
i
i
i天需要的人.
我们用志愿者代表流量,我们要用志愿者覆盖
[
l
,
r
+
1
]
[l,r+1]
[l,r+1]之间的连边,因为覆盖
(
i
,
i
+
1
)
(i,i+1)
(i,i+1)的边通过志愿者流过了至少
a
i
a_i
ai,此时我们把
(
i
,
i
+
1
)
(i,i+1)
(i,i+1)的容量上限定义为
−
a
i
-a_i
−ai,然后求总流量为0的最小费用流即可.由于网络流容量非负,所以加上
I
N
F
INF
INF即可.
#include <bits/stdc++.h>
#define fi first
#define se second
#define lc (x << 1)
#define rc (x << 1 | 1)
#define gc getchar() //(p1==p2&&(p2=(p1=buf)+fread(buf,1,size,stdin),p1==p2)?EOF:*p1++)
#define mk make_pair
#define pii pair<int, int>
#define pll pair<ll, ll>
#define pb push_back
#define IT iterator
#define V vector
#define TP template <class o>
#define TPP template <typename t1, typename t2>
#define SZ(a) ((int)a.size())
#define all(a) a.begin(), a.end()
#define rep(i, a, b) for (auto i = a; i <= b; i++)
#define REP(i, a, b) for (auto i = b; i >= a; i--)
#define debug(x) cerr << #x << ' ' << '=' << ' ' << x << endl
using namespace std;
typedef unsigned ui;
typedef long long ll;
typedef unsigned long long ull;
typedef double db;
typedef long double ld;
const int N = 10010, size = 1 << 20, mod = 998244353, inf = 0x7fffffff;
const ll INF = 1e15;
// char buf[size],*p1=buf,*p2=buf;
TP void qr(o& x) {
char c = gc;
x = 0;
int f = 1;
while (!isdigit(c)) {
if (c == '-')
f = -1;
c = gc;
}
while (isdigit(c))
x = x * 10 + c - '0', c = gc;
x *= f;
}
TP void qw(o x) {
if (x / 10)
qw(x / 10);
putchar(x % 10 + '0');
}
TP void pr1(o x) {
if (x < 0)
x = -x, putchar('-');
qw(x);
putchar(' ');
}
TP void pr2(o x) {
if (x < 0)
x = -x, putchar('-');
qw(x);
putchar(10);
}
// math
ll gcd(ll a, ll b) { return !a ? b : gcd(b % a, a); }
ll lcm(ll a, ll b) { return a / gcd(a, b) * b; }
ll power(ll a, ll b = mod - 2, ll p = mod) {
ll c = 1;
while (b) {
if (b & 1)
c = c * a % p;
b /= 2;
a = a * a % p;
}
return c;
}
TP void cmax(o& x, o y) {
if (x < y)
x = y;
}
void cmax(int& x, int y) { x = x - y >> 31 ? y : x; }
TP void cmin(o& x, o y) {
if (x > y)
x = y;
}
void cmin(int& x, int y) { x = x - y >> 31 ? x : y; }
TPP void ad(t1& x, t2 y) {
x += y;
if (x >= mod)
x -= mod;
}
TPP void dl(t1& x, t2 y) {
x -= y;
if (x < 0)
x += mod;
}
int n, m, a[N];
namespace FLOW {
struct edge {
int y, next, c, d;
} a[N * 4];
int len = 1, last[N];
void ins(int x, int y, int c, int d) {
a[++len] = (edge){ y, last[x], c, d };
last[x] = len;
}
void add(int x, int y, int c, int d) {
ins(x, y, c, d);
ins(y, x, 0, -d);
}
int q[N], d[N], st, ed, pre[N], f[N];
bool vis[N];
void EK() {
int ans = 0;
while (1) {
int l = 1, r = 2;
q[1] = st;
f[st] = inf;
f[ed] = 0;
memset(d, 63, sizeof d);
d[st] = 0;
while (l ^ r) {
int x = q[l++];
vis[x] = 0;
if (l == N)
l = 1;
for (int k = last[x], y; k; k = a[k].next) {
y = a[k].y;
if (a[k].d + d[x] < d[y] && a[k].c) {
f[y] = min(f[x], a[k].c);
d[y] = d[x] + a[k].d;
pre[y] = k;
if (!vis[y]) {
vis[y] = 1;
q[r++] = y;
if (r == N)
r = 1;
}
}
}
}
if (!f[ed]) break;
ans += f[ed] * d[ed];
int x = ed, k, t = f[ed];
while (x ^ st) {
k = pre[x];
a[k].c -= t;
a[k ^ 1].c += t;
x = a[k ^ 1].y;
}
}
pr2(ans);
}
}
using namespace FLOW;
void solve() {
qr(n);
qr(m);
st = 0;
ed = n + 2;
add(st, 1, inf, 0);
add(n + 1, ed, inf, 0);
int l, r, x;
rep(i, 1, n) qr(x), add(i, i + 1, inf - x, 0);
while (m--) {
qr(l);
qr(r);
qr(x);
add(l, ++r, inf, x);
}
EK();
}
int main() {
int T = 1;
// qr(T);
while (T--)
solve();
return 0;
}
线性规划的套路转化网络流方法
先写成松弛型,设有
m
m
m个约束,
n
n
n个随机变量.
我们在前后加入两个约束
0
=
0
0=0
0=0.
然后,我们进行初等行变换/转对偶形式.(本题直接差分即可)
最后每个约束形如 若干变量+常数 = 0.
且对于每个变量而言,只会在约束中出现两次,并且系数为
±
1
\pm 1
±1.
我们可以把每一个约束看作一个点,那么约束本身可以视作流量守恒.
我们把每一个随机变量看作一条边的流量(
x
i
≥
0
x_i\ge 0
xi≥0),然后对于
−
1
-1
−1看作出边,
+
1
+1
+1看作入边.
例如,
x
i
x_i
xi存在于
j
,
k
j,k
j,k两行,系数分别为
−
1
,
+
1
-1,+1
−1,+1,那么就相当于给连一条边
(
j
,
k
)
(j,k)
(j,k).
在本题中由于
x
x
x带权,我们直接给
x
x
x相关的边带上权跑最小费用最大流即可.
code
#3550. [ONTAK2010]Vacation#3550. [ONTAK2010]Vacation
令
x
i
=
0
o
r
1
x_i=0~or ~1
xi=0 or 1,表示第
i
i
i个值选还是不选.
m
a
x
i
m
i
z
e
∑
i
=
1
n
x
i
a
i
maximize \sum_{i=1}^n x_ia_i
maximize∑i=1nxiai
s
.
t
.
∑
i
=
l
l
+
n
−
1
x
i
≤
k
,
x
i
≤
1
,
x
i
≥
0
s.t. \sum_{i=l}^{l+n-1} x_i\le k,x_i\le 1,x_i\ge 0
s.t.∑i=ll+n−1xi≤k,xi≤1,xi≥0.
容易看出这是一个全幺模矩阵…
然后直接套板子…
const int N = 610, M = 1010;
const db eps = 1e-8, inf = 1e10;
int n, m, type, id[N + M], tp[N];
db a[M][N];
void pivot(int r, int c) {
swap(id[r + n], id[c]);
db t = -a[r][c];
a[r][c] = -1;
rep(j, 0, n) a[r][j] /= t;
rep(i, 0, m) if (i ^ r && a[i][c]) {
t = a[i][c];
a[i][c] = 0;
rep(j, 0, n) a[i][j] += a[r][j] * t;
}
}
void solve() {
FOR(i, n) id[i] = i;
db t;
while (1) {
int i = 0, j = 0;
db w = eps;
FOR(k, n) if (a[0][k] > w) w = a[0][j = k];
if (!j)
break;
w = inf;
FOR(k, m) if (-a[k][j] > eps && (t = -a[k][0] / a[k][j]) < w) w = t, i = k;
if (!i) {
puts("Unbounded");
return;
}
pivot(i, j);
}
pr2((int)(a[0][0] + 0.2)), a[0][0] = 0;
rep(i, n + 1, n + m) tp[id[i]] = i - n;
if (type)
FOR(i, n) printf("%.9lf ", a[tp[i]][0]);
}
int main() {
qr(n);
int k;
qr(k);
FOR(i, 3 * n) scanf("%lf", &a[0][i]);
FOR(i, 2 * n + 1) {
a[++m][0] = k;
rep(j, i, i + n - 1) a[m][j] = -1;
}
n *= 3;
FOR(i, n) a[++m][0] = 1, a[m][i] = -1;
solve();
return 0;
}
当然我们还是可以转为网络流跑的…
写出所有约束:
{
0
=
0
∑
i
=
l
l
+
n
−
1
x
i
+
D
l
=
k
(
l
∈
[
1
,
2
n
+
1
]
)
0
=
0
\begin{cases} 0=0\\ \sum_{i=l}^{l+n-1} x_i+D_l=k(l\in [1,2n+1])\\ 0=0 \end{cases}
⎩⎪⎨⎪⎧0=0∑i=ll+n−1xi+Dl=k(l∈[1,2n+1])0=0
差分后得到:
{
x
1
+
x
2
+
.
.
.
+
x
n
+
D
1
−
k
=
0
x
i
+
n
−
1
−
x
i
−
1
+
D
i
−
D
i
−
1
=
0
(
i
∈
[
2
,
2
n
+
1
]
)
k
−
∑
i
=
2
n
+
1
3
n
x
i
−
D
2
n
+
1
=
0
\begin{cases} x_1+x_2+...+x_n+D_1-k=0\\ x_{i+n-1}-x_{i-1}+D_i-D_{i-1}=0(i\in [2,2n+1])\\ k-\sum_{i=2n+1}^{3n} x_i-D_{2n+1}=0 \end{cases}
⎩⎪⎨⎪⎧x1+x2+...+xn+D1−k=0xi+n−1−xi−1+Di−Di−1=0(i∈[2,2n+1])k−∑i=2n+13nxi−D2n+1=0
设 ( x , y , c , d ) (x,y,c,d) (x,y,c,d)表示弧 ( x , y ) (x,y) (x,y)的容量为c,费用为d.
则我们可以这样连边.
i
∈
[
2
n
+
1
,
3
n
]
,
x
i
:
(
2
n
+
2
,
i
−
n
+
1
,
1
,
a
i
)
i\in [2n+1,3n] ,x_i:(2n+2,i-n+1,1,a_i)
i∈[2n+1,3n],xi:(2n+2,i−n+1,1,ai)注意到
x
i
x_i
xi的上限为1,所以我们把它设为容量上限.
i
∈
[
1
,
2
n
]
,
x
i
:
(
i
+
1
,
max
(
1
,
i
−
n
+
1
)
,
1
,
a
i
)
i\in [1,2n],x_i:(i+1,\max(1,i-n+1),1,a_i)
i∈[1,2n],xi:(i+1,max(1,i−n+1),1,ai).
i
∈
[
1
,
2
n
+
1
]
,
D
i
:
(
i
+
1
,
i
,
I
N
F
,
0
)
i\in[1,2n+1],D_i:(i+1,i,INF,0)
i∈[1,2n+1],Di:(i+1,i,INF,0).
(
s
t
,
2
n
+
2
,
k
,
0
)
,
(
1
,
e
d
,
k
,
0
)
(st,2n+2,k,0),(1,ed,k,0)
(st,2n+2,k,0),(1,ed,k,0).
qr(n);
int k;
qr(k);
st = 2*n+3;
ed = st+1;
add(st, 2 * n + 2, k, 0);
add(1, ed, k, 0);
rep(i, 1, 3 * n) qr(m), add(min(2 * n + 2, i+1), max(i - n + 1, 1), 1, m);
rep(i, 1, 2 * n + 1) add(i+1, i, inf, 0);
EK();
P3337 [ZJOI2013]防守战线
比前面的更套路…
const int N = 1010, M = 10010;
const db eps = 1e-8, inf = 1e10;
int n, m;
db a[M][N];
void pivot(int r, int c) {
db t = -a[r][c];
a[r][c] = -1;
rep(j, 0, n) a[r][j] /= t;
rep(i, 0, m) if (i ^ r && a[i][c]) {
t = a[i][c];
a[i][c] = 0;
rep(j, 0, n) a[i][j] += t * a[r][j];
}
}
void solve() {
db t;
while (1) {
int i = 0, j = 0;
db w = -eps;
FOR(k, m) if (a[k][0] < w) w = a[i = k][0];
if (!i)
break;
FOR(k, n) if (a[i][k] > eps) {
j = k;
break;
}
pivot(i, j);
}
while (1) {
int i = 0, j = 0;
db w = eps;
FOR(k, n) if (a[0][k] > w) w = a[0][j = k];
if (!j)
break;
w = inf;
FOR(k, m) if (-a[k][j] > eps && (t = -a[k][0] / a[k][j]) < w) w = t, i = k;
pivot(i, j);
}
pr2((int)(-a[0][0] + 0.2));
}
int main() {
qr(n);
qr(m);
FOR(i, n) qr(a[0][i]), a[0][i] *= -1;
FOR(i, m) {
int l, r, c;
qr(l);
qr(r);
qr(c);
a[i][0] = -c;
rep(j, l, r) a[i][j] = 1;
}
solve();
return 0;
}
没事转化成对偶线性规划也行~
我们还是尝试用线性规划转网络流(这样的练习是防止比赛的时候卡单纯形的空间…
设
s
i
s_i
si为前
i
i
i个位置的总塔数.
线性规划为:
{
m
i
n
i
m
i
z
e
∑
i
=
1
n
d
i
(
s
i
−
s
i
−
1
)
=
∑
i
=
1
n
s
i
(
d
i
−
d
i
+
1
)
s
.
t
.
s
i
−
s
i
−
1
≥
0
,
i
∈
[
1
,
n
]
s
r
j
−
s
l
j
−
1
≥
c
j
\begin{cases} minimize &\sum_{i=1}^n d_i(s_i-s_{i-1})=\sum_{i=1}^n s_i(d_i-d_{i+1})\\ s.t. &s_i-s_{i-1}\ge 0,i\in[1,n]\\ & s_{r_j}-s_{l_j-1}\ge c_j\\ \end{cases}
⎩⎪⎨⎪⎧minimizes.t.∑i=1ndi(si−si−1)=∑i=1nsi(di−di+1)si−si−1≥0,i∈[1,n]srj−slj−1≥cj
我们必须写成对偶形式,因为这样每一个随机变量均出现两次,才能视作流量.
对偶线性规划:
令
p
i
p_i
pi对应前
n
+
1
n+1
n+1个约束(转对偶之前加入
0
≥
0
0\ge 0
0≥0),
q
i
q_i
qi对应后
m
m
m约束.
{
m
a
x
m
i
z
e
∑
j
=
1
m
c
j
q
j
p
i
−
p
i
+
1
+
∑
r
j
=
i
q
j
−
∑
l
j
−
1
=
i
q
j
≤
d
i
−
d
i
+
1
\begin{cases} maxmize & \sum_{j=1}^m c_jq_j\\ &p_i-p_{i+1} +\sum_{r_j=i} q_j-\sum_{l_j-1=i} q_j\le d_i-d_{i+1} \\ \end{cases}
{maxmize∑j=1mcjqjpi−pi+1+∑rj=iqj−∑lj−1=iqj≤di−di+1
容易发现对式子累加得到
0
≤
0
0\le 0
0≤0,所以所有条件都必须取等.
连边:
p
i
:
(
i
,
i
+
1
,
I
N
F
,
0
)
p_i:(i,i+1,INF,0)
pi:(i,i+1,INF,0).
q
j
:
(
l
j
−
1
,
r
j
,
I
N
F
,
c
j
)
q_j:(l_j-1,r_j,INF,c_j)
qj:(lj−1,rj,INF,cj).
d
i
+
1
>
d
i
,
(
s
t
,
i
,
d
i
+
1
−
d
i
,
0
)
,
e
l
s
e
(
i
,
e
d
,
d
i
−
d
i
+
1
,
0
)
d_{i+1}>d_i,(st,i,d_{i+1}-d_i,0),else ~(i,ed,d_i-d_{i+1},0)
di+1>di,(st,i,di+1−di,0),else (i,ed,di−di+1,0)
4363
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
