bzoj3571/洛谷P3236/loj2205 画框 KM算法

最新推荐文章于 2022-08-18 11:30:02 发布

litble

最新推荐文章于 2022-08-18 11:30:02 发布

阅读量368

点赞数 1

分类专栏：计算几何图论文章标签： KM算法叉积最小乘积生成树

本文链接：https://blog.csdn.net/litble/article/details/79679391

版权

图论同时被 2 个专栏收录

16 篇文章 0 订阅

订阅专栏

计算几何

6 篇文章 0 订阅

订阅专栏

题目分析

真TM神题。

对于一种匹配方案，我们将其记为一个点 $\sum A_{i,p_i} ,\sum B_{i,p_i})$ ，那么我们要求横纵坐标相乘最小的一个点。乘积相等的两点，一定在同一条反比例函数曲线上。反比例函数曲线，绝对值越小，越靠近坐标轴，所以我们不停维护答案的话，应该要维护出一个下凸包。

考虑这么一种算法，首先找到x坐标最小的点A和y坐标最小的点B，然后每次找一个在直线AB左下方的最远点C，在递归处理AC左下方和CB左下方，直到找不到为止。

那么如何找这个C呢？发现由于AB的长度是定值，那么三角形ABC的面积又可以反应C到直线AB的距离大小。而向量AC和AB的叉积又等于这个三角形面积的2倍，所以我们要求最大叉积。（当然，如果C在直线AB上方，那么叉积就是个负数）

$AC \times AB=(X_C-X_A)(Y_B-Y_A)-(X_B-X_A)(Y_C-Y_A)=$
$Y_B-Y_A)X_C+(X_A-X_B)Y_C-X_AY_B+X_BY_A$

后面的定值暂且不管，看向前面，我们只要构造一个二分图，一边是画一边是画框，其中i与j之间的边权是 $Y_B-Y_A)a_{i,j}+(X_A-X_B)b_{i,j}$ ，然后做一个带权二分图匹配即可。

但是这题丧心病狂卡费用流，所以用KM算法吧。

~~复杂度O(跑得过)~~

代码

#include<bits/stdc++.h>
using namespace std;
#define RI register int
int read() {
	int q=0;char ch=' ';
	while(ch<'0'||ch>'9') ch=getchar();
	while(ch>='0'&&ch<='9') q=q*10+ch-'0',ch=getchar();
	return q;
}
const int N=71,inf=0x3f3f3f3f;
struct point{int x,y;};
int a[N][N],b[N][N],g[N][N];
int visa[N],visb[N],cp[N],dl[N],expa[N],expb[N];
int T,n;

int dfs(int x) {
	visa[x]=1;
	for(RI i=1;i<=n;++i) {
		if(visb[i]) continue;
		int kl=expa[x]+expb[i]-g[x][i];
		if(kl==0) {
			visb[i]=1;
			if(!cp[i]||dfs(cp[i])) {cp[i]=x;return 1;}
		}
		else dl[i]=min(dl[i],kl);
	}
	return 0;
}
point km() {
	for(RI i=1;i<=n;++i) expb[i]=cp[i]=0;
	for(RI i=1;i<=n;++i) {
		expa[i]=-inf;
		for(RI j=1;j<=n;++j) expa[i]=max(expa[i],g[i][j]);
	}
	for(RI i=1;i<=n;++i) {
		 for(RI j=1;j<=n;++j) dl[j]=inf;
		 while("niconiconi") {
		 	for(RI j=1;j<=n;++j) visa[j]=visb[j]=0;
		 	if(dfs(i)) break;
		 	int kl=inf;
		 	for(RI j=1;j<=n;++j) if(!visb[j]) kl=min(kl,dl[j]);
		 	for(RI j=1;j<=n;++j) {
		 		if(visa[j]) expa[j]-=kl;
		 		if(visb[j]) expb[j]+=kl;
		 		else dl[j]-=kl;
		 	}
		 }
	}
	point re=(point){0,0};
	for(RI i=1;i<=n;++i) re.x+=a[cp[i]][i],re.y+=b[cp[i]][i];
	return re;
}
int work(point l,point r) {
	for(RI i=1;i<=n;++i)
		for(RI j=1;j<=n;++j) g[i][j]=a[i][j]*(r.y-l.y)+b[i][j]*(l.x-r.x);
	point mid=km();
	if((l.x==mid.x&&l.y==mid.y)||(r.x==mid.x&&r.y==mid.y))
		return min(l.x*l.y,r.x*r.y);
	else return min(work(l,mid),work(mid,r));
}
int main()
{
	T=read();
	while(T--) {
		n=read();
		for(RI i=1;i<=n;++i)
			for(RI j=1;j<=n;++j) a[i][j]=read();
		for(RI i=1;i<=n;++i)
			for(RI j=1;j<=n;++j) b[i][j]=read();
		for(RI i=1;i<=n;++i)
			for(RI j=1;j<=n;++j) g[i][j]=-a[i][j];
		point L=km();
		for(RI i=1;i<=n;++i)
			for(RI j=1;j<=n;++j) g[i][j]=-b[i][j];
		point R=km();
		printf("%d\n",work(L,R));
	}
    return 0;
}