AT_geocon2013_b 玉座の間 题解

前置知识：费用流

题目分析

题意

平面上有 \(n\) 个点，你要给每个点移动位置，使得最终的图形沿 \(y\) 轴对称。求移动的欧几里得距离之和的最小值。

其实从题干想到费用流确实很神奇。

但是想到之后后继步骤还是很简单的。

解法

首先易得以下几条结论：

• 把一个点移到 \(y\) 轴的另一面并不优。（感性理解一下即可）

• 把一个点移到 \(y\) 轴上也是满足条件的。

所以可能会想到分成 \(y\) 轴左侧和 \(y\) 轴右侧两部分进行匹配。

但是两边都可以互相匹配。

所以得到下一个神仙想法：拆点。

将第 \(i\) 个点拆为 \(L_i\) 和 \(R_i\) 两个点。

源点 \(S\) 向 \(L_i\) 连流量为 \(1\)，费用为 \(0\) 的边。

\(R_i\) 向 汇点 \(T\) 连流量为 \(1\)，费用为 \(0\) 的边。

然后对于第 \(i\) 个点来说，可以枚举另一个点 \(j\)，然后分成三种情况连边：

• \(i=j\)，自己和自己匹配，相当于移到 \(y\) 轴上。此时从 \(L_i\) 向 \(R_i\) 连流量为 \(1\)，费用为 \(x_i\) 的边。

• \(i\) 和 \(j\) 在 \(y\) 轴同一侧，即 \(x_i\times x_j\leq0\) 时。因为不优，直接跳过。

• \(i\) 与 \(j\) 不在 \(y\) 轴同一侧，即 \(x_i\times x_j>0\) 时。此时相当于将自己挪到 \(j\) 沿 \(y\) 轴对称的点上。从 \(L_i\) 向 \(R_j\) 连流量为 \(1\)，费用为 \(\cfrac{\sqrt{(x_i+x_j)^2+(y_i-y_j)^2}}{2}\) 的边。

最后一种情况费用为 \(\cfrac{\sqrt{(x_i+x_j)^2+(y_i-y_j)^2}}{2}\) 而非 \(\sqrt{(x_i+x_j)^2+(y_i-y_j)^2}\)，是因为匹配是两边同时匹配，意味着 \(i\) 匹配 \(j\)，同时 \(j\) 也会匹配 \(i\)。匹配两次，所以费用是距离的 \(\cfrac{1}{2}\)。

最后跑最小费用最大流，结束。

Code

一开始浮点数的费用流出了点问题。

后面找到了，判费用时用的是 !cost[i]

所以我把费用乘以 \(10^8\)，计算完后再除以 \(10^8\) 即可。（只保留 7 位小数）

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

template<typename Tp, size_t sizn, size_t sizm>
struct costflow
{
    #define pb __gnu_pbds
    typedef pb::priority_queue<pair<Tp, int>, 
                               greater<pair<Tp, int> >, 
                               pb::pairing_heap_tag> pairing_heap;

    int cnt=1, s=sizn-2, t=sizn-3;
    Tp delta=0, MXflow=0, MNcost=0;
    void link(int u, int v, Tp w, Tp f) 
    {
        to [++cnt]=v;       cap [cnt]=w;
        nxt[ cnt ]=head[u]; head[ u ]=cnt;
        cost[cnt ]=f;       from[cnt]=u;
        to [++cnt]=u;       cap [cnt]=0;
        nxt[ cnt ]=head[v]; head[ v ]=cnt;
        cost[cnt ]=-f;      from[cnt]=v;
    }

    private:

    bitset<sizn> vis;
    pairing_heap pq; 
    typename pairing_heap :: point_iterator it[sizn];
    Tp val[sizm<<1], dis[sizn], flow[sizn], cost[sizm<<1], cap[sizm<<1];
    int head[sizn], to[sizm<<1], nxt[sizm<<1], now[sizm<<1], from[sizm<<1];

    bool SPFA()
    {
        queue<int> q;
        memset(dis, 0x3f, sizeof dis);
        vis.reset();
        dis[t]=0; q.push(t);
        vis[t]=1;
        while (!q.empty())
        {
            int u=q.front(); q.pop();
            vis[u]=0;
            for(int i=head[u];i;i=nxt[i])
            {
                int v=to[i];
                Tp f=val[i^1], c=cap[i^1], len=cost[i^1];
                if(f<c&&dis[v]>dis[u]+len)
                {
                    dis[v]=dis[u]+len;
                    if(!vis[v]) vis[v]=1, q.push(v);
                }
            }
        }
        return dis[s]!=dis[0];
    }

    void reduce()
    {
        for(int i=2;i<=cnt;i++) cost[i]+=dis[to[i]]-dis[from[i]];
        delta+=dis[s];
    }

    bool Dijkstra() 
    {
        memset(dis, 0x3f, sizeof dis);
        memset(it, 0, sizeof it);
        dis[t] = 0; 
        it[t] = pq.push({dis[t], t});
        while(!pq.empty()) 
        {
            int u = pq.top().second; pq.pop();
            for(int j=head[u];j;j=nxt[j]) 
            {
                int v=to[j];
                Tp f=val[j^1], c=cap[j^1], len=cost[j^1];
                if(f<c&&dis[v]>dis[u]+len) 
                {
                    dis[v]=dis[u]+len;
                    if(it[v]==NULL) it[v]=pq.push({dis[v], v});
                    else pq.modify(it[v], {dis[v], v});
                }
            }
        }
        return dis[s] != dis[0];
    }
    
    // dfs 找增广路并处理
    Tp dfs(int idx, Tp sum)
    {
        if(idx==t) return sum;
        vis[idx]=1;
        Tp k, ret=sum;
        for(int i=head[idx];i&&sum;i=nxt[i])
        {
            int v=to[i];
            Tp f=val[i], c=cap[i], len=cost[i];
            if(!vis[v]&&f<c&&!len)
            {
                k=dfs(v, min(ret, c-f));
                val[i]+=k;
                val[i^1]-=k;    
                ret-=k;
            }
        }
        return sum-ret;
    }

    void augment()
    {
        Tp curflow=0;
        vis.reset();
        while((curflow=dfs(s, dis[0])))
        {
            MXflow+=curflow;
            MNcost+=curflow*delta;
            vis.reset();
        }
    }
    public:
    void PrimalDual()
    {
        if(!SPFA()) return;
        reduce(); augment();
        while(Dijkstra()) {reduce(); augment();}
    }
};

costflow<long long, 10005, 100005> cf;

struct dot
{
    int x, y;
    dot(int X=0, int Y=0) : x(X), y(Y) {}
    double operator-(const dot &b) { return sqrt((x-b.x)*(x-b.x)+(y-b.y)*(y-b.y)); }
    dot rev() { return {-x, y}; }
};

dot vtc[105];

#define l(x) x<<1
#define r(x) x<<1|1

int main()
{
    int n;
    cin>>n;
    for(int i=1;i<=n;i++)
        cin>>vtc[i].x>>vtc[i].y;
    for(int i=1;i<=n;i++)
    {
        cf.link(cf.s, l(i), 1, 0);
        cf.link(r(i), cf.t, 1, 0);
        cf.link(l(i), r(i), 1, abs(vtc[i].x)*1e8);
        for(int j=1;j<=n;j++)
        {
            if(vtc[j].x*vtc[i].x>=0) continue;
            cf.link(l(i), r(j), 1, (vtc[i]-(vtc[j].rev()))/2*1e8);
        }
    }
    cf.PrimalDual();
    printf("%.7f", (double)cf.MNcost/1e8);
}