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\le 0$ 时。因为不优，直接跳过。

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

最后一种情况费用为 $\frac{\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$。匹配两次，所以费用是距离的 $\frac{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);
}