bzoj3571

同样的最小乘积XXX模型，这里显然是二分图带权匹配

我不会写KM……于是写了个费用流，由于是稠密图，会退化到n^4

然后本地跑了56s，交上去过了………………一定是我电脑太慢……

改天写个KM吧

const inf=14000*14000;
type node=record
po,next,flow,cost:longint;
end;
point=record
x,y:longint;
end;

var e:array[0..2000010] of node;
ma,pre,cur,p,d:array[0..145] of longint;
v:array[0..145] of boolean;
a,b,c:array[0..75,0..75] of longint;
q:array[0..2000010] of longint;
j,len,i,tt,t,n,ans:longint;
p1,p2:point;

procedure add(x,y,f,c:longint);
begin
inc(len);
e[len].po:=y;
e[len].flow:=f;
e[len].cost:=c;
e[len].next:=p[x];
p[x]:=len;
end;

procedure build(x,y,f,c:longint);
begin
add(x,y,f,c);
add(y,x,0,-c);
end;

function spfa:boolean;
var f,r,x,y,i,j:longint;
begin
fillchar(v,sizeof(v),false);
for i:=1 to t do
d[i]:=inf;
d[0]:=0;
f:=1;
r:=1;
q[1]:=0;
while f<=r do
begin
x:=q[f];
v[x]:=false;
i:=p[x];
while i<>-1 do
begin
y:=e[i].po;
if e[i].flow>0 then
if d[x]+e[i].cost<d[y] then
begin
d[y]:=d[x]+e[i].cost;
cur[y]:=i;
pre[y]:=x;
if not v[y] then
begin
inc(r);
q[r]:=y;
v[y]:=true;
end;
end;
i:=e[i].next;
end;
inc(f);
end;
if d[t]=inf then exit(false) else exit(true);
end;

procedure change(j:longint);
begin
dec(e[j].flow);
inc(e[j xor 1].flow);
end;

function get:point;
var i,j,x,y:longint;
begin
len:=-1;
fillchar(p,sizeof(p),255);
for i:=1 to n do
begin
build(0,i,1,0);
build(i+n,t,1,0);
end;
x:=1; y:=1;
for i:=1 to n do
for j:=1 to n do
begin
build(i,j+n,1,c[i,j]);
if c[i,j]<c[x,y] then
begin
x:=i;
y:=j;
end;
end;
ma[y+n]:=x;
i:=p[x];
while i<>-1 do
begin
if (e[i].po=y+n) then change(i);
if e[i].po=0 then change(i xor 1);
i:=e[i].next;
end;
i:=p[y+n];
while i<>-1 do
begin
if e[i].po=t then
begin
change(i);
break;
end;
i:=e[i].next;
end;
while spfa do
begin
i:=t;
while i<>0 do
begin
j:=cur[i];
change(j);
if (i>n) and (i<>t) then
ma[i]:=pre[i];
i:=pre[i];
end;
end;
get.x:=0; get.y:=0;
for i:=1 to n do
begin
j:=ma[i+n];
// write(j,' ');
inc(get.x,a[j,i]);
inc(get.y,b[j,i]);
end;
// writeln;
if get.x*get.y<ans then ans:=get.x*get.y;
end;

function cross(a,b,c:point):longint;
begin
exit((a.x-c.x)*(b.y-c.y)-(a.y-c.y)*(b.x-c.x));
end;

procedure work(p1,p2:point);
var p:point;
i,j:longint;
begin
for i:=1 to n do
for j:=1 to n do
c[i,j]:=a[i,j]*(p1.y-p2.y)+b[i,j]*(p2.x-p1.x);
p:=get;
if cross(p2,p,p1)>=0 then exit;
work(p1,p);
work(p,p2);
end;

begin
readln(tt);
while tt>0 do
begin
dec(tt);
readln(n);
t:=2*n+1;
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
c[i,j]:=a[i,j];
ans:=14000*140000;
p1:=get;
for i:=1 to n do
for j:=1 to n do
c[i,j]:=b[i,j];
p2:=get;
work(p1,p2);
writeln(ans);
end;
end.

View Code