[BZOJ3571][HNOI2014]画框-最小乘积KM

画框

Description

小T准备在家里摆放几幅画，为此他买来了N幅画和N个画框。为了体现他的品味，小T希望能合理地搭配画与画框，使得其显得既不过于平庸也不太违和。对于第 幅画与第 个画框的配对，小T都给出了这个配对的平凡度Aij 与违和度Bij 。整个搭配方案的总体不和谐度为每对画与画框平凡度之和与每对画与画框违和度的乘积。具体来说，设搭配方案中第i幅画与第Pi个画框配对，则总体不和谐度为

小T希望知道通过搭配能得到的最小的总体不和谐度是多少。

Input

输入文件第 行是一个正整数T ，表示数据组数，接下来是T组数据。
对于每组数据，第 行是一个正整数N，表示有N对画和画框。
第2到第N+1行，每行有N个非负整数，第i+1 行第j个数表示Aij 。
第N+2到第2*N+1行，每行有N个非负整数，第i+N+1 行第j个数表示Bij 。

Output

包含T行，每行一个整数，表示最小的总体不和谐度

Sample Input

1
3
4 3 2
2 3 4
3 2 1
2 3 2
2 2 4
1 1 3

Sample Output

30

HINT

第1幅画搭配第3个画框，第2幅画搭配第1个画框，第3 幅画搭配第2个画框，则总体不和谐度为30

N<=70,T<=3,Aij<=200,Bij<=200

---

最小乘积生成树的进化体……
表示还需要学习一个

---

思路:
考虑使用最小乘积生成树的套路:

对于一种方案，令$x=\sum a$，$y= \sum b$。
那么，考虑使用点对$(x,y)$来描述一种匹配方案，此时每个方案的答案为$x*y$。
可以发现，把所有方案在平面上列出，则可能的最优解构成一个下凸壳。
那么求出这个凸壳即可~

问题在于，点数过多，无法像通常一样全部求出再去进行凸包构造。

此时可以考虑另一种构造方案:
首先令边权为$a_{ij}$，跑一边最小权匹配，得到点$A$，然后令边权为$b_{ij}$再跑一遍，得到点$B$。
可以发现，这两个点构成了凸壳的边界，因为实际最优解一定在这两个节点与原点的区域内。

然后，每次递归地去寻找距离上一层传下来的两个点在原点方向上最远的一个节点$C$，再将$(A,C)$和$(C,B)$分别继续向下递归，直到向量$AB$和向量$AC$的叉积大于$0$再停下。
最后就能得到一个凸壳了~

考虑如何找到这个$C$。
由于$C$在原点方向上距离$AB$最远，那么有$S_{\triangle}ABC$最大。
于是由于$AB$、$AC$两个向量的叉积的面积的一半等同于$S_{\triangle}ABC$，那么就是要最大化
$(B_x-A_x)*(C_y-A_y)-(B_y-A_y)*(C_x-A_x)$
拆一发:
$(B.x-A.x)C_y+(A_y-B_y)C_x-(B.x-A.x)A_y+(B_y-A_y)A_x$
由于不与$C$相关的部分为常数，因此目标是最大化前半部分。

那么令图上每一条边的权值为$(B_x-A_x)*b_{ij}+(A_y-B_y)*a_{ij}$，跑一遍最小权匹配即可得到$C$!

最小权匹配可以考虑使用边权取负后的$KM$算法实现~

#include<bits/stdc++.h>
using namespace std;

inline int read()
{
    int x=0;char ch=getchar();
    while(ch<'0' || '9'<ch)ch=getchar();
    while('0'<=ch && ch<='9')x=x*10+(ch^48),ch=getchar();
    return x;
}

inline void chkmin(int &a,int b){if(a>b)a=b;}
inline void chkmax(int &a,int b){if(a<b)a=b;}

struct pr
{
    int x,y;
    pr(int _x=0,int _y=0){x=_x;y=_y;}
    inline pr operator - (pr b)
    {
        return pr(x-b.x,y-b.y);
    }
};
inline int cross(pr a,pr b){return a.x*b.y-a.y*b.x;}

const int N=79;
int n,ans;
int g[N][N],a[N][N],b[N][N];
int lx[N],ly[N],match[N],rec[N];
bool visx[N],visy[N];

inline bool dfs(int x)
{
    if(visx[x])return 0;
    visx[x]=1;
    for(int i=1;i<=n;i++)
        if(!visy[i] && g[x][i]==lx[x]+ly[i])
        {
            visy[i]=1;
            if(!match[i] || dfs(match[i]))
                return match[i]=x,1;
        }
        else if(!visy[i])
            chkmin(rec[i],lx[x]+ly[i]-g[x][i]);
    return 0;
}

inline pr km()
{
    memset(ly,0,sizeof(ly[0])*(n+9));
    memset(lx,128,sizeof(lx[0])*(n+9));
    memset(match,0,sizeof(match[0])*(n+9));
    for(int i=1;i<=n;i++)
        for(int j=1;j<=n;j++)
            chkmax(lx[i],g[i][j]);
    for(int v=1;v<=n;v++)
    {
        memset(visx,0,sizeof(bool)*(n+9));
        memset(visy,0,sizeof(bool)*(n+9));
        memset(rec,127,sizeof(int)*(n+9));
        while(!dfs(v))
        {
            int mv=1e9+7;
            for(int i=1;i<=n;i++)
                if(!visy[i])
                    chkmin(mv,rec[i]);
            for(int i=1;i<=n;i++)
            {
                if(visx[i])lx[i]-=mv,visx[i]=0;
                if(visy[i])ly[i]+=mv,visy[i]=0;
            }
        }
    }

    int reta=0,retb=0;
    for(int i=1;i<=n;i++)
    {
        reta+=a[match[i]][i];
        retb+=b[match[i]][i];
    }

    chkmin(ans,reta*retb);
    return pr(reta,retb);
}

inline void work(pr pa,pr pb)
{
    for(int i=1;i<=n;i++)
        for(int j=1;j<=n;j++)
            g[i][j]=-((pb.x-pa.x)*b[i][j]+(pa.y-pb.y)*a[i][j]);
    pr c=km();
    if(cross(pb-pa,c-pa)>=0)return;
    work(pa,c);work(c,pb);
}

int mina()
{
    n=read();ans=1e9+7;
    for(int i=1;i<=n;i++)
        for(int j=1;j<=n;j++)
            a[i][j]=read();
    for(int i=1;i<=n;i++)
        for(int j=1;j<=n;j++)
            b[i][j]=read();

    for(int i=1;i<=n;i++)
        for(int j=1;j<=n;j++)
            g[i][j]=-a[i][j];
    pr a=km();
    for(int i=1;i<=n;i++)
        for(int j=1;j<=n;j++)
            g[i][j]=-b[i][j];
    pr b=km();

    work(a,b);
    printf("%d\n",ans);
    return 0;
}

int main()
{
    for(int T=read();T;T--)
        mina();

    return 0;
}