BZOJ 3571

传送门

首先,$N$个配$N$个还带权值,那就显然是用KM了

这道题的思路和最小乘积生成树的思路差不多,可惜BZOJ上一道裸最小乘积生成树被权限了…

考虑把每一种方案抽象成一个点$(\sum A,\sum B)$,然后把所有点扔到一个坐标系里面,然后答案就是使得一个反比例函数$x\cdot y=k$的$k$最小的点

拿KM找出$\sum A$最小的点和$\sum B$最小的点(带权最小匹配就把边权取反算带权最大匹配),最优解一定优于这两个点(或位于两个点之一),也就是在这两个点连成的线段下方

设最优解为$c$,$\sum A$最小的点为$a$,$\sum B$最小的点为$b$

易证明$c$一定在下凸壳上

那么题目就成了求解$a$到$b$的下凸壳。那么问题来了,总点数太多,不可能得到所有的点,怎么求下凸壳?

有一个显然的结论,$\vec{ab}$下方距离$\vec{ab}$最远的点一定在下凸壳上,即使得$S\Delta abc$最大的$c$一定在下凸壳上。

$$S\Delta abc=\frac{\vec{ab}\times\vec{ac}}{2}$$
展开再合起来能得到 (向量箭头不一样高差评)
$$S\Delta abc=\frac{\vec{ab}\times\vec{c}+\vec{a}\times\vec{b}}{2}$$

后面那坨已知,扔掉,即求$\vec{ab}\times\vec{c}$的最大值

找到$c$之后,更新答案,$c$不一定是最优解,我们要将下凸壳遍历,因此继续左右分治下去。

分析时间复杂度,由于要将下凸壳全遍历,设下凸壳上的点数为$S$,每次${\rm KM}$的复杂度是$\Theta(N^3)$,总复杂度是$\Theta(SN^3)$。

在实践中,这个算法是很优秀的,但是还是存在很坏的情况……(跑的还没有暴力快)


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
program bzoj_3571;
const inf=maxlongint>>1;
type point=record x,y:longint; end;
var lx,ly,slack,pair:array[1..70]of longint;
vx,vy:array[1..70]of boolean;
w,a,b:array[1..70,1..70]of longint;
t,n,i,j:longint;
pa,pb:point;
function min(a,b:longint):longint;
begin
if a<b then exit(a) else exit(b);
end;
function find(x:longint):boolean;
var t,y:longint;
begin
vx[x]:=true;
for y:=1 to n do
begin
if vy[y] then continue;
t:=lx[x]+ly[y]-w[x,y];
if t=0 then
begin
vy[y]:=true;
if (pair[y]=0)or(find(pair[y])) then
begin
pair[y]:=x;
exit(true);
end;
end else if t<slack[y] then slack[y]:=t;
end;
exit(false);
end;
function KM:point;
var i,j,d:longint;
begin
fillchar(pair,sizeof(pair),0);
for i:=1 to n do lx[i]:=-inf;
fillchar(ly,sizeof(ly),0);
for i:=1 to n do
for j:=1 to n do if w[i,j]>lx[i] then lx[i]:=w[i,j];
for i:=1 to n do
begin
fillchar(slack,sizeof(slack),$3f);
repeat
fillchar(vx,sizeof(vx),false);
fillchar(vy,sizeof(vy),false);
if find(i) then break;
d:=inf;
for j:=1 to n do if (not vy[j])and(slack[j]<d) then d:=slack[j];
for j:=1 to n do if vx[j] then dec(lx[j],d);
for j:=1 to n do if vy[j] then inc(ly[j],d) else dec(slack[j],d);
until false;
end;
KM.x:=0;
KM.y:=0;
for i:=1 to n do inc(KM.x,a[pair[i],i]);
for i:=1 to n do inc(KM.y,b[pair[i],i]);
end;
function solve(l,r:point):longint;
var t:point;
i,j:longint;
begin
for i:=1 to n do
for j:=1 to n do w[i,j]:=a[i,j]*(r.y-l.y)-b[i,j]*(r.x-l.x);
t:=KM;
if ((t.x=l.x)and(t.y=l.y))or((t.x=r.x)and(t.y=r.y)) then exit(min(l.x*l.y,r.x*r.y));
exit(min(solve(l,t),solve(t,r)));
end;
begin
readln(t);
for t:=1 to t do
begin
readln(n);
for i:=1 to n do
for j:=1 to n do read(a[i,j]);
for i:=1 to n do
for j:=1 to n do read(b[i,j]);
for i:=1 to n do
for j:=1 to n do w[i,j]:=-a[i,j];
pa:=KM;
for i:=1 to n do
for j:=1 to n do w[i,j]:=-b[i,j];
pb:=KM;
writeln(solve(pa,pb));
end;
end.