线性规划与单纯形法

标准型

  • 最大化 (sumlimits_{j=1}^{n}c_jx_j)
  • 满足约束 (sumlimits_{j=1}^{n}a_{ij}x_jle b_i) (i=1,2,...,m)
    • ​ (x_jge 0) (j=1,2,...,m)

松弛型

  • 最大化 (sumlimits_{j=1}^{n}c_jx_j)
  • 满足约束 (x_{i+n}=b_i-sumlimits_{j=1}^{n}a_{ij}x_j)
    • ​ (x_jge 0) (j=1,2,...,n+m)
  • 线性规划的解空间是一个凸性区域，局部最优解就是全局最优解。

单纯形

基变量:在松弛型等式左侧的所有变量。

非基变量:在松弛型等式右侧的所有变量。

  • 如果常数项都非负，我们可以拿出一个基本解，即所有基变量的值为右侧常数项，非基变量的值为(0)。

转轴(pivot)操作

  • 选择一个基变量(x_B)和一个非基变量(x_N)，将其互换(称这个非基变量为换入变量，基变量为换出变量)。
  • (x_B=b_i-sumlimits_{j=1}^{n}a_{ij}x_j)
  • (x_N=(b_i-sumlimits_{j eq N}a_{ij}x_j-x_B)/a_{iN})
  • 将目标函数以及其余所有限制中(x_N)用这个等式替代。

simplex 操作

  • 从基本解出发，每次找一个目标函数系数为正的非基变量和一个对他约束最强的基变量进行转轴操作。
  • 保证每次转轴操作都会使得目标函数增大。

initialization操作

  • 初始化，即找到一组基本解。

  • 做法:引入一个辅助线性规划:

    • 最大化 (-x_0)
    • 满足约束 (x_{i+n}=b_i-sumlimits_{j=1}^{n}a_{ij}x_j+x_0)
      • ​ (x_jge 0)

    如果原线性规划存在一个解为((x_1,...,x_{n+m}))，那么((0,x_1,..,x_{n+m}))为辅助线性规划的一个可行解。 我们最大化(-x_0)那么最优解就是这个可行解。

  • 构造辅助线性规划的初始解，我们把(x_0)当做换入变量，找到最小的(b_l)，把(x_{l+n})当做换出变量执行一次转轴操作。

    • (x_0=-b_l+sumlimits_{j=1}^{n}a_{lj}x_j+x{l+n})

  • 对于其他的约束，变为

    • (x_{i+n}=-b_l+b_i+sumlimits_{j=1}(a_{lj}-a{ij})x_j+x{l+n})
    • 操作之后辅助线性规划有一个可行的初始基本解了，执行(simplex)操作即可。

时间复杂度

(O(simplex) imes O(nm)=O(???)=O(可过))

可以证明单纯形法一定能够在有限步之后终止，但是最坏情况算法的时间复杂度为指数级别的。

对偶问题

给定一个原始线性规划:

  • 最小化 (c^{T}x)
  • 满足约束 (Ax ge b)
    (x ge 0)

其对偶线性规划为:

  • 最大化 (b^{T}y)

  • 满足约束 (A^{T}y le c)
    (y ge 0)

线性规划弱对偶性:若(x=(x_1,...x_n))是原问题的一个可行解，(y=(y1,...ym))是对偶问题的可行解。那么:

  • (sumlimits_{j=1}^{n}c_jx_jge sumlimits_{i=1}^{m}b_i y_i)

线性规划对偶性:若(x=(x_1,...x_n))是原问题的一个最优解，(y=(y1,...ym))是对偶问题的最优解。那么:

  • (sumlimits_{j=1}^{n}c_jx_j=sumlimits_{i=1}^{m}b_iy_i)

根据这个性质可以将一系列问题的初始化操作去掉，转化成对偶后的线性规划进行求解。

网络流解线性规划问题

https://www.cnblogs.com/suika/p/9079315.html

  • 对于一些线性规划问题，我们通常能够转化成 每个变量的都出现两次，且系数分别为(+1)和(-1)。
  • 就是这样的模型，我们可以用网络流的方法巧妙的解决。
  • 首先网络流有性质：对于除了源点和汇点的其他点，有 流入的总量等于流出的总量。
  • 这让我们有一个思路，把点当成每一个方程，把边看成给予限制的变量，对于每个变量，把(+1)看成流出(-1)看成流入，找到对应方程的位置连边。

常见的问题转化成线性规划

最短路:

  • 最大化 (d_t)
  • 满足约束 (d_s=0 d_{vi}le d_{ui}+l_i)

最大流

  • 设(f_{u,v})表示(u)到(v)的实际流量。
  • 最大化 (sumlimits_{u}f_{s,u})
  • 满足约束 (f_{u,v}=f_{v,u})
    • ​ (sumlimits_{i eq s,t}f_{u,i}=0)

全幺模矩阵

  • 全幺模矩阵是指任何一个行数和列数相等的满秩子矩阵行列式的值都是(1)或(-1)的矩阵。
  • 对于一个线性规划中的(A)为全幺模矩阵，那么所有约束系数都是(1,0,-1)中的一个。
  • 也就是说，直接执行线性规划就能得到相同的整数规划问题的正确答案。

模板

http://uoj.ac/problem/179

#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;
#define N 105
#define eps 1e-8
typedef double f2;
f2 a[N][N];
int n,m,type,id[N],tp[N];
void pivot(int r,int c) {
	swap(id[r+n],id[c]);
	f2 tmp=-a[r][c];
	a[r][c]=-1;
	int i,j;
	for(i=0;i<=n;i++) a[r][i]/=tmp;
	for(i=0;i<=m;i++) if(a[i][c]&&i!=r) {
		tmp=a[i][c]; a[i][c]=0;
		for(j=0;j<=n;j++) a[i][j]+=tmp*a[r][j];
	}
}
void solve() {
	int i,j,k;
	for(i=1;i<=n;i++) id[i]=i;
	while(1) {
		i=j=0; 
		f2 tmp=-eps;
		for(k=1;k<=m;k++) if(a[k][0]<tmp) tmp=a[i=k][0];
		if(!i) break;
		for(k=1;k<=n;k++) if(a[i][k]>eps) {j=k; break;}
		if(!j) {puts("Infeasible"); return ;}
		pivot(i,j);
	}
	while(1) {
		i=j=0; 
		f2 tmp=eps;
		for(k=1;k<=n;k++) if(a[0][k]>tmp) tmp=a[0][j=k];
		if(!j) break;
		tmp=1e9;
		for(k=1;k<=m;k++) if(a[k][j]<-eps&&-a[k][0]/a[k][j]<tmp) {
			i=k; tmp=-a[k][0]/a[k][j];
		}
		if(!i) {puts("Unbounded"); return ;}
		pivot(i,j);
	}
	printf("%.9f
",a[0][0]);
	for(i=n+1;i<=n+m;i++) tp[id[i]]=i-n;
	if(type) {
		for(i=1;i<=n;i++) printf("%.9f ",tp[i]?a[tp[i]][0]:0);
	}
}
int main() {
	scanf("%d%d%d",&n,&m,&type);
	int i,j;
	for(i=1;i<=n;i++) scanf("%lf",&a[0][i]);
	for(i=1;i<=m;i++) {
		for(j=1;j<=n;j++) scanf("%lf",&a[i][j]),a[i][j]*=-1; 
		scanf("%lf",&a[i][0]);
	}
	solve();
}

论文中的例题

https://lydsy.com/JudgeOnline/problem.php?id=3265 (fake)


https://community.topcoder.com/stat?c=problem_statement&pm=8605&rd=12012


BZOJ1061: [Noi2008]志愿者招募

https://lydsy.com/JudgeOnline/problem.php?id=1061

题解: https://www.cnblogs.com/suika/p/9079315.html

代码:

#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;
#define N 10050
#define M 100050
#define S (10049)
#define T (10048)
#define inf (1<<30)
int head[N],to[M],nxt[M],flow[M],val[M],cnt=1,path[N],Q[N],l,r,dis[N],inq[N],n,m,a[N],c[M],s[N],t[N];
inline void add(int u,int v,int f,int w) {
    to[++cnt]=v; nxt[cnt]=head[u]; head[u]=cnt; flow[cnt]=f; val[cnt]=w;
    to[++cnt]=u; nxt[cnt]=head[v]; head[v]=cnt; flow[cnt]=0; val[cnt]=-w;
}
bool spfa() {
    memset(dis,0x3f,sizeof(dis));
    memset(path,0,sizeof(path));
    l=r=0; Q[r++]=S; dis[S]=0; inq[S]=1;
    while(l!=r) {
        int x=Q[l++],i; if(l==n+10) l=0; inq[x]=0;
        for(i=head[x];i;i=nxt[i]) {
            if(flow[i]&&dis[to[i]]>dis[x]+val[i]) {
                dis[to[i]]=dis[x]+val[i];
                path[to[i]]=i^1;
                if(!inq[to[i]]) {
                    Q[r++]=to[i]; if(r==n+10) r=0; inq[to[i]]=1;
                }
            }
        }
    }
    return path[T];
}
void mcmf() {
    int minc=0,maxf=0;
    while(spfa()) {
        int nf=1<<30,i;
        for(i=T;i!=S;i=to[path[i]]) {
            nf=min(nf,flow[path[i]^1]);
        }
        for(i=T;i!=S;i=to[path[i]]) {
            flow[path[i]]+=nf;
            flow[path[i]^1]-=nf;
            minc+=nf*val[path[i]^1];
        }
        maxf+=nf;
    }
    printf("%d
",minc);
}
int main() {
    scanf("%d%d",&n,&m);
    int i;
    for(i=1;i<=n;i++) {
        scanf("%d",&c[i]); add(i+1,i,inf,0);
    }
    for(i=1;i<=m;i++) {
        scanf("%d%d%d",&s[i],&t[i],&a[i]);
        add(s[i],t[i]+1,inf,a[i]);
    }
    for(i=1;i<=n+1;i++) {
        if(c[i]-c[i-1]>0) {
            add(S,i,c[i]-c[i-1],0);
        }else {
            add(i,T,c[i-1]-c[i],0);
        }
    }
    mcmf();
}

BZOJ3550[ONTAK2010]Vacation&&BZOJ1283:序列

https://lydsy.com/JudgeOnline/problem.php?id=3550

题解:https://www.cnblogs.com/suika/p/9080922.html

代码:

#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;
#define N 1050
#define M 100050
#define S (n-m+3)
#define T (n-m+4)
#define inf (1<<30)
int head[N],to[M],nxt[M],cnt=1,path[N],dis[N],Q[N],l,r,inq[N],flow[M],val[M];
int n,m,K,a[N];
inline void add(int u,int v,int f,int w) {
    to[++cnt]=v; nxt[cnt]=head[u]; head[u]=cnt; flow[cnt]=f; val[cnt]=w;
    to[++cnt]=u; nxt[cnt]=head[v]; head[v]=cnt; flow[cnt]=0; val[cnt]=-w;
}
bool spfa() {
    memset(dis,0x80,sizeof(dis));
    memset(path,0,sizeof(path));
    l=r=0; Q[r++]=S; dis[S]=0; inq[S]=1;
    while(l!=r) {
        int x=Q[l++],i; if(l==T+1) l=0; inq[x]=0;
        for(i=head[x];i;i=nxt[i]) {
            if(flow[i]&&dis[to[i]]<dis[x]+val[i]) {
                dis[to[i]]=dis[x]+val[i];
                path[to[i]]=i^1;
                if(!inq[to[i]]) {
                    inq[to[i]]=1; Q[r++]=to[i]; if(r==T+1) r=0;
                }
            }
        }
    }
    return path[T];
}
void mcmf() {
    int minc=0,maxf=0,i;
    while(spfa()) {
        int nf=1<<30;
        for(i=T;i!=S;i=to[path[i]]) {
            nf=min(nf,flow[path[i]^1]);
        }
        for(i=T;i!=S;i=to[path[i]]) {
            flow[path[i]]+=nf;
            flow[path[i]^1]-=nf;
            minc+=nf*val[path[i]^1];
        }
        maxf+=nf;
    }
    printf("%d
",minc);
}
int main() {
    scanf("%d%d%d",&n,&m,&K);
    int i;
    for(i=1;i<=n;i++) {
        scanf("%d",&a[i]);
    }
    add(S,1,K,0); add(n-m+2,T,K,0);
    for(i=1;i<=n-m+1;i++) add(i,i+1,inf,0);
    for(i=1;i<=n;i++) {
        if(i<=m) add(1,i+1,1,a[i]);
        else if(i>n-m) add(i-m+1,n-m+2,1,a[i]);
        else add(i-m+1,i+1,1,a[i]);
    }
    mcmf();
}

BZOJ1937: [Shoi2004]Mst 最小生成树

https://lydsy.com/JudgeOnline/problem.php?id=1937

同BZOJ3118。

空间开成([9900][770])可过。

需要判负(0)的情况。

代码:

#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
using namespace std;
typedef double f2;
#define N 55
#define M 770
//maxm=766
#define eps 1e-8
int head[N],to[N<<1],nxt[N<<1],cnt,n,m,fa[N],val[N<<1];
int idx[N],dep[N],tot;
f2 a[9950][M];
struct A {
	int x,y,k,z,ist;
	void read() {
		scanf("%d%d%d",&x,&y,&z);
	}
}e[M];
inline void add(int u,int v,int w) {
	to[++cnt]=v; nxt[cnt]=head[u]; head[u]=cnt; val[cnt]=w;
}
void dfs(int x,int y) {
	int i;
	fa[x]=y;
	dep[x]=dep[y]+1;
	for(i=head[x];i;i=nxt[i]) if(to[i]!=y) {
		dfs(to[i],x);
		idx[to[i]]=val[i];
	}
}
int lca(int x,int y) {
	while(x!=y) {
		if(dep[x]<dep[y]) y=fa[y];
		else x=fa[x];
	}
	return x;
}
void insert(int t,int nt) {
	tot++;
	a[tot][t]=1;
	a[tot][nt]=1;
	a[tot][0]=e[nt].z-e[t].z;
}
void build(int id) {
	int x=e[id].x,y=e[id].y;
	int l=lca(x,y);
	while(x!=l) {
		insert(idx[x],id);
		x=fa[x];
	}
	while(y!=l) {
		insert(idx[y],id);
		y=fa[y];
	}
}
void pivot(int r,int c,int n,int m) {
	f2 t=-a[r][c]; 
	a[r][c]=-1;
	int i,j;
	for(i=0;i<=n;i++) a[r][i]/=t;
	for(i=0;i<=m;i++) if(a[i][c]&&i!=r) {
		t=a[i][c];
		a[i][c]=0;
		for(j=0;j<=n;j++) a[i][j]+=t*a[r][j];
	}
}
void simplex(int n,int m) {
	int i,j,k;
	f2 t;
	while(1) {
		i=j=0;
		t=-eps;
		for(k=1;k<=m;k++) if(a[k][0]<t) t=a[i=k][0];
		if(!i) break;
		for(k=1;k<=n;k++) if(a[i][k]>eps) {j=k; break;}
		// if(!j) {puts("NO"); return ;}
		pivot(i,j,n,m);
	}
	
	while(1) {
		i=j=0;
		t=eps;
		for(k=1;k<=n;k++) if(a[0][k]>t) t=a[0][j=k];
		if(!j) break;
		t=1e9;
		for(k=1;k<=m;k++) if(a[k][j]<-eps&&-a[k][0]/a[k][j]<t) {
			i=k, t=-a[k][0]/a[k][j];
		}
		// if(!i) {puts("INF"); return ;}
		pivot(i,j,n,m);
	}
	if(fabs(a[0][0])<eps) puts("0");
	else printf("%.0f
",-a[0][0]);
}
int ID[N][N];
int main() {
	scanf("%d%d",&n,&m);
	int i;
	for(i=1;i<=m;i++) {
		e[i].read();
		if(e[i].x>e[i].y) swap(e[i].x,e[i].y);
		ID[e[i].x][e[i].y]=i;
	}
	int x,y;
	for(i=1;i<n;i++) {
		scanf("%d%d",&x,&y);
		if(x>y) swap(x,y);
		e[ID[x][y]].ist=1;
	}
	for(i=1;i<=m;i++) {
		if(e[i].ist) {
			add(e[i].x, e[i].y, i);
			add(e[i].y, e[i].x, i);
		}
	}
	dfs(1,0);
	for(i=1;i<=m;i++) {
		if(e[i].ist==0) {
			build(i);
			a[0][i]=-1;
		}else {
			a[0][i]=-1;
		}
	}
	simplex(m,tot);
}

https://www.codechef.com/FEB12/problems/FLYDIST


http://poj.org/problem?id=3689


https://www.codechef.com/JUNE15/problems/CHEFBOOK


BZOJ3118: Orz the MST

https://lydsy.com/JudgeOnline/problem.php?id=3118

分析:

  • 考虑一条非树边会约束哪些树边。
  • 然后就可以做了，注意求最小要把(c)都取负数求最大。

代码:

#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;
typedef double f2;
#define N 305
#define M 1005
#define eps 1e-8
int head[N],to[N<<1],nxt[N<<1],cnt,n,m,fa[N],val[N<<1];
int idx[N],dep[N],tot;
f2 a[4050][M];
struct A {
	int x,y,k,z,a,b,ist;
	void read() {
		scanf("%d%d%d%d%d%d",&x,&y,&z,&ist,&a,&b);
	}
}e[M];
inline void add(int u,int v,int w) {
	to[++cnt]=v; nxt[cnt]=head[u]; head[u]=cnt; val[cnt]=w;
}
void dfs(int x,int y) {
	int i;
	fa[x]=y;
	dep[x]=dep[y]+1;
	for(i=head[x];i;i=nxt[i]) if(to[i]!=y) {
		dfs(to[i],x);
		idx[to[i]]=val[i];
	}
}
int lca(int x,int y) {
	while(x!=y) {
		if(dep[x]<dep[y]) y=fa[y];
		else x=fa[x];
	}
	return x;
}
void insert(int t,int nt) {
	tot++;
	a[tot][t]=1;
	a[tot][nt]=1;
	a[tot][0]=e[nt].z-e[t].z;
}
void build(int id) {
	int x=e[id].x,y=e[id].y;
	int l=lca(x,y);
	while(x!=l) {
		insert(idx[x],id);
		x=fa[x];
	}
	while(y!=l) {
		insert(idx[y],id);
		y=fa[y];
	}
}
void pivot(int r,int c,int n,int m) {
	f2 t=-a[r][c]; 
	a[r][c]=-1;
	int i,j;
	for(i=0;i<=n;i++) a[r][i]/=t;
	for(i=0;i<=m;i++) if(a[i][c]&&i!=r) {
		t=a[i][c];
		a[i][c]=0;
		for(j=0;j<=n;j++) a[i][j]+=t*a[r][j];
	}
}
void simplex(int n,int m) {
	int i,j,k;
	f2 t;
	while(1) {
		i=j=0;
		t=-eps;
		for(k=1;k<=m;k++) if(a[k][0]<t) t=a[i=k][0];
		if(!i) break;
		for(k=1;k<=n;k++) if(a[i][k]>eps) {j=k; break;}
		// if(!j) {puts("NO"); return ;}
		pivot(i,j,n,m);
	}
	
	while(1) {
		i=j=0;
		t=eps;
		for(k=1;k<=n;k++) if(a[0][k]>t) t=a[0][j=k];
		if(!j) break;
		t=1e9;
		for(k=1;k<=m;k++) if(a[k][j]<-eps&&-a[k][0]/a[k][j]<t) {
			i=k, t=-a[k][0]/a[k][j];
		}
		// if(!i) {puts("INF"); return ;}
		pivot(i,j,n,m);
	}
	printf("%.0f
",-a[0][0]);
}
int main() {
	scanf("%d%d",&n,&m);
	int i;
	for(i=1;i<=m;i++) {
		e[i].read();
		if(e[i].ist) {
			add(e[i].x, e[i].y, i);
			add(e[i].y, e[i].x, i);
		}
	}
	dfs(1,0);
	for(i=1;i<=m;i++) {
		if(e[i].ist==0) {
			build(i);
			a[0][i]=-e[i].a;
		}else {
			a[0][i]=-e[i].b;
		}
	}
	simplex(m,tot);
}

BZOJ3112: [Zjoi2013]防守战线

https://lydsy.com/JudgeOnline/problem.php?id=3112

这道题和BZOJ1061是同一类题。

我这里为了练单纯形法而写了这个。

可以看出单纯形法效率之高。

代码:

#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
using namespace std;
typedef double f2;
typedef long double f3;
#define N 1050
#define M 10050
#define eps 1e-9
#define inf 1e10
f2 a[M][N];
int n,m,x;
void pivot(int r,int c,int n,int m) {
	f2 t=-a[r][c];
	a[r][c]=-1;
	int i,j;
	for(i=0;i<=n;i++) a[r][i]/=t;
	for(i=0;i<=m;i++) if(abs(a[i][c])>eps&&i!=r) {
		t=a[i][c]; a[i][c]=0;
		for(j=0;j<=n;j++) a[i][j]+=t*a[r][j];
	}
}
void simplex(int n,int m) {
	int i,j,k;
	f2 t;
	while(1) {
		i=j=0;
		t=-eps;
		for(k=1;k<=m;k++) if(a[k][0]<t) t=a[i=k][0];
		if(!i) break;
		for(k=1;k<=n;k++) if(a[i][k]>eps) {j=k; break;}
		if(!j) {puts("NO"); return ;}
		pivot(i,j,n,m);
	}
	while(1) {
		i=j=0;
		t=eps;
		for(k=1;k<=n;k++) if(a[0][k]>t) t=a[0][j=k];
		if(!j) break;
		t=inf;
		for(k=1;k<=m;k++) if(a[k][j]<-eps&&(-a[k][0]/a[k][j]<t)) {
			t=-a[k][0]/a[k][j]; i=k;
		}
		if(!i) {puts("INF"); return ;}
		pivot(i,j,n,m);
	}
}
int main() {
	scanf("%d%d",&n,&m);
	int i,x,l,r,j;
	for(i=1;i<=n;i++) {
		scanf("%d",&x);
		a[0][i]=-x;
	}
	for(i=1;i<=m;i++) {
		scanf("%d%d",&l,&r);
		for(j=l;j<=r;j++) {
			a[i][j]=1;
		}
		scanf("%lf",&a[i][0]);
		a[i][0]*=-1;
	}
	simplex(n,m);
	if(abs(a[0][0])<eps) puts("0");
	else printf("%.0f
",-a[0][0]);
}

https://vjudge.net/problem/UVA-10498


http://codeforces.com/contest/375/problem/E

我的做法是设(i)节点(x_i)表示是否被交换过。
那么如果原来是(0)，现在就是(x_i)，原来是1，现在就是(1-x_i)。
需要注意(x_ile 1)。

代码:

#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
using namespace std;
typedef long long ll;
typedef double f2;
#define N 515
#define eps 1e-8
int head[N],to[N<<1],nxt[N<<1],val[N<<1],n,cnt,wa[N];
ll dis[N][N],X;
f2 a[N<<1][N];
inline void add(int u,int v,int w) {
	to[++cnt]=v; nxt[cnt]=head[u]; head[u]=cnt; val[cnt]=w;
}
void dfs(int x,int y,int rt,ll d) {
	dis[rt][x]=d;
	int i;
	for(i=head[x];i;i=nxt[i]) if(to[i]!=y) {
		dfs(to[i],x,rt,d+val[i]);
	}
}
void pivot(int r,int c,int n,int m) {
	int i,j;
	f2 t=-a[r][c]; a[r][c]=-1;
	for(i=0;i<=n;i++) a[r][i]/=t;
	for(i=0;i<=m;i++) if(abs(a[i][c])>eps&&i!=r) {
		t=a[i][c]; a[i][c]=0;
		for(j=0;j<=n;j++) a[i][j]+=t*a[r][j];
	}
}
void simplex(int n,int m) {
	int i,j,k;
	f2 t;
	while(1) {
		i=j=0;
		t=-eps;
		for(k=1;k<=m;k++) if(a[k][0]<t) t=a[i=k][0];
		if(!i) break;
		for(k=1;k<=n;k++) if(a[i][k]>eps) {
			j=k; break;
		}
		if(!j) {puts("-1"); return ;}
		pivot(i,j,n,m);
	}
	while(1) {
		i=j=0;
		t=eps;
		for(k=1;k<=n;k++) if(a[0][k]>t) t=a[0][j=k];
		if(!j) break;
		t=1e9;
		for(k=1;k<=m;k++) if(a[k][j]<-eps&&(-a[k][0]/a[k][j])<t) {
			t=-a[k][0]/a[k][j]; i=k;
		}
		if(!i) {puts("INF"); return ;}
		pivot(i,j,n,m);
	}
	if(abs(a[0][0])<eps) puts("0");
	else if(-a[0][0]>n) puts("-1");
	else printf("%.0f
",-a[0][0]/2);
}
int main() {
	scanf("%d%lld",&n,&X);
	int i,x,y,z,j,sum=0;
	for(i=1;i<=n;i++) scanf("%d",&wa[i]),sum+=wa[i];
	for(i=1;i<n;i++) {
		scanf("%d%d%d",&x,&y,&z);
		add(x,y,z); add(y,x,z);
	}
	for(i=1;i<=n;i++) {
		dfs(i,0,i,0);
	}
	for(i=1;i<=n;i++) a[0][i]=-1;
	for(i=1;i<=n;i++) {
		int t=0;
		for(j=1;j<=n;j++) {
			if(dis[i][j]<=X) {
				if(wa[j]) t++,a[i][j]=-1;
				else a[i][j]=1;
			}
		}
		a[i][0]=t-1;
	}
	int t=0;
	for(i=1;i<=n;i++) {
		if(wa[i]) t++,a[n+1][i]=-1;
		else a[n+1][i]=1;
	}
	a[n+1][0]=t-sum;
	t=0;
	for(i=1;i<=n;i++) {
		if(wa[i]) t--,a[n+2][i]=1;
		else a[n+2][i]=-1;
	}
	a[n+2][0]=t+sum;
	int tot=n+2;
	for(i=1;i<=n;i++) {
		a[++tot][i]=-1;
		a[tot][0]=1;
	}
	simplex(n,tot);
}

https://community.topcoder.com/stat?c=problem_statement&pm=14049&rd=16652

【推广】 免费学中医，健康全家人
原文地址：https://www.cnblogs.com/suika/p/10022545.html