【算法】最小乘积生成树 & 最小乘积匹配 （HNOI2014画框）

今天考试的时候果然题目太难于是我就放弃了……转而学习了一下最小乘积生成树。

最小乘积生成树定义：

（摘自网上一篇博文）。

我们主要解决的问题就是当k = 2时，如何获得最小的权值乘积。我们注意到一张图可以有很多棵生成树，我们将每一棵生成树的权值记为（x, y），表示第一种权值之和为x, 第二种权值之和为y. 这样，很自然联想到二维平面上的坐标，每一棵生成树即为这个平面上的一个点。我们所想要寻找的点就是x * y最小的点。这样的点在什么位置？显然，若x1 <= x2, y1 <= y2,1号点的权值必然更小。所以我们的答案只可能处于这张平面图上的凸包的下凸壳上。

于是我们找到A，B两点，一个离y轴最近，一个离x轴最近，这两个点一定是下凸壳的两个端点。之后，我们再寻找到与AB距离最远的点C，用点C 更新答案后再以AC，BC为新的边向下递归求解。此时问题来了：如何找到这一个距离最远，且在AB下方的C点呢？我们将距离转化为面积，使用叉积求解。因为要求C点在AB下方，所以得到的叉积必为负数。又因为|叉积| = 四边形面积，所以得到的叉积必然是负的值中绝对值最大的那一个，即求解出与AB构成的叉积最小的C点。

然后就开始考虑式子的转化：min (B - A) * (C - A) = (B.x - A.x) (C.x - A.x) - (C.x - A.x) (B.y - A.y); 化开这个式子，省去常数部分，我们发现所求就是（A.y - B.y）* a[i][j] - （A.x - B.x）* b[i][j] 最小。我们考虑将这个东西看做权值，就可以用Kruskal求出使这个值最小的C点了。如果是匹配的话，则将i --> j 视作匹配的权值，将权值取反（因为要求求最小）后跑KM算法获得最大权值匹配。

下面的代码是仿照着的，但觉得写的很漂亮，放在这里大家可以参考一下。感谢原本的博主~

#include <bits/stdc++.h>
using namespace std;
#define maxn 1000
#define INF 99999999
int n, ans = INF, lx[maxn], ly[maxn], s[maxn], match[maxn];
int T, g[maxn][maxn], a[maxn][maxn], b[maxn][maxn];
bool visx[maxn], visy[maxn];

struct vec
{
    int x, y;
};

vec operator -(vec a, vec b)
{
    return (vec) { b.x - a.x, b.y - a.y };
}

int operator *(vec a, vec b)
{
    return a.x * b.y - a.y * b.x;
}

int read()
{
    int x = 0, k = 1;
    char c;
    c = getchar();
    while(c < '0' || c > '9') { if(c == '-') k = -1; c = getchar(); }
    while(c >= '0' && c <= '9') x = x * 10 + c - '0', c = getchar();
    return x * k;
}

struct Graph
{
    void build(int wx, int wy)
    {
        for(int i = 1; i <= n; i ++)
            for(int j = 1; j <= n; j ++)
                g[i][j] = - (wx * a[i][j] + wy * b[i][j]);
    }    
    
    bool dfs(int u)
    {
        visx[u] = 1;
        for(int v = 1; v <= n; v ++)
        {
            if(visy[v]) continue;
            int tem = lx[u] + ly[v] - g[u][v];
            if(!tem)
            {
                visy[v] = 1;
                if(!match[v] || dfs(match[v]))
                {
                    match[v] = u;
                    return 1;
                }
            }
            else s[v] = min(s[v], tem);
        }
        return false;
    }
    
    vec KM()
    {
        memset(lx, 0, sizeof(lx)), memset(ly, 0, sizeof(ly));
        memset(match, 0, sizeof(match));
        for(int i = 1; i <= n; i ++)
            for(int j = 1; j <= n; j ++)
                lx[i] = max(lx[i], g[i][j]);
        for(int i = 1; i <= n; i ++)
        {
            memset(s, 63, sizeof(s));
            while(23333)
            {
                memset(visx, 0, sizeof(visx)), memset(visy, 0, sizeof(visy));
                if(dfs(i)) break;
                int tem = INF;
                for(int j = 1; j <= n; j ++)
                    if(!visy[j]) tem = min(tem, s[j]);
                for(int j = 1; j <= n; j ++)
                    if(visx[j]) lx[j] -= tem;
                for(int j = 1; j <= n; j ++)
                    if(visy[j]) ly[j] += tem;
                    else s[j] -= tem;
            }
        }
        vec re; re.x = 0, re.y = 0;
        for(int i = 1; i <= n; i ++)
            re.x += a[match[i]][i], re.y += b[match[i]][i];
        return re;
    }
}G;

void Solve(vec A, vec B)
{
    G.build(A.y - B.y, B.x - A.x);
    vec C = G.KM();
    ans = min(ans, C.x * C.y);
    if((A - B) * (A - C) >= 0) return;
    Solve(A, C), Solve(C, B);
}

int main()
{
    T = read();
    while(T --)
    {
        n = read();
        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();
        G.build(1, 0); 
        vec A = G.KM();
        G.build(0, 1);
        vec B = G.KM();
        ans = min(A.x * A.y, B.x * B.y);
        Solve(A, B);
        printf("%d\n", ans);
    }
    return 0;
}