【网络流】最大密度子图

参考：https://wenku.baidu.com/view/986baf00b52acfc789ebc9a9.html

目录

简介
原理
代码

简介

给定无向图 \(G=(V,E)\) ，其子图记为 \(G'=(V',E')\) ，在所有子图构成的集合中，密度 \(D=\frac{|E'|}{|V'|}\) 最大的元素称为最大密度子图。

原理

如何解决这样的问题呢？

首先我们进行二分答案，二分出一个 \(g\) ，对于 \(\max(|E'|-g|V'|)\) ，如果其大于 \(0\) ，那么说明答案更大，否则说明答案更小，如此下去便可求出答案。

故下面我们来考察 \(\max(|E'|-g|V'|)\) ，方便起见，将它写成一个以 \(g\) 为自变量的函数： \(h(g) = |E'|-g|V'|\) ，最大化 \(h(g)\) 即可。

注意到目前的式子没有太好的性质（可以和网络流中的最小割建立联系），考虑进一步变换。

对于 \(|E'|\) ，我们有：

\[|E'| = \frac{\sum_{u\in V'}d_u - cnt[V',\overline{V'}]}{2} \]

其中 \(d_u\) 表示 \(u\) 的度数，\(cnt[V',\overline{V'}]\) 表示非子图 \(G'\) 内部边的个数，也就是 \(cnt[V',\overline{V'}]\) 之间相连的边数。

我们有 \(h(g) = \frac{\sum_{u\in V'}d_u - cnt[V',\overline{V'}]}{2}-g|V'|\) ，整理一下

\[h(g) = \frac{1}{2}(\sum_{u\in V'}d_u - cnt[V',\overline{V'}]-2g|V'|) \]

而 \(2g|V'| = \sum_{u\in V'}2g\) ，进一步有：

\[h(g) = \frac{1}{2}(\sum_{u\in V'}{(d_u-2g)} - cnt[V',\overline{V'}]) \]

因为需要和最小割建立联系，所以问题变成最小化 \(-h(g)\) ：

\[-h(g) = \frac{1}{2}(\sum_{u\in V'}{(2g-d_u)} + cnt[V',\overline{V'}]) \]

根据式子的特征，我们这样建图：\(V\) 中的点之间连接容量为 \(1\) 的边，\(V\) 中的点 \(u\) 向汇点 \(t\) 连接容量为 \(2g-d_u+U\) 的边（\(U\) 为一个足够保证容量非负的数），源点 \(s\) 向 \(V\) 中的点 \(u\) 连接容量为 \(U\) 的边。

我们将从接下来的公式推导中看出如此建图是合理而且正确的。

注意到割的容量由且只由 \(s \rightarrow \overline{V'}\) 、\(V' \rightarrow t\) 以及 \(V'\rightarrow \overline{V'}\) 构成。

因此对于一个割，其容量满足：

\[c[s,t]=\sum_{u\in \overline{V'}}U + \sum_{u\in V'}(2g-d_u+U) + \sum_{u\in V'}\sum_{v\in \overline{V'}}c(u,v) \]

进一步化为：

\[c[s,t]=|V|U + \sum_{u\in V'}(2g-d_u) + \sum_{u\in V'}\sum_{v\in \overline{V'}}c(u,v) \]

而 \(\sum_{u\in V'}(-d_u) + \sum_{u\in V'}\sum_{v\in \overline{V'}}c(u,v)\) 恰好为 \(-2|E'|\) ，故上式可进而得到：

\[c[s,t]=|V|U + \sum_{u\in V'}(2g) - 2|E'| \]

整理一下，有：

\[c[s,t]=|V|U + 2(g|V'| - |E'|) \]

这意味着最小割对应着最小的 \(-h(g) = g|V'|-|E'|\) ，也就是最大的 \(h(g) = |E'|-g|V'|\) 。

至此，问题解决。

那么在具体实现的时候， \(U\) 应该如何选取呢？

注意到 \(U\) 是为了保证 \(2g-d_u\) 非负，所以取 \(U = |E'|\) 即可。

模板题传送门：https://www.acwing.com/activity/content/problem/content/2624/1/

代码

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

const int N=110, M=1000+2*N<<1;
const double INF=1e10, eps=1e-8;

struct edge{
    int u, v;
}edges[M];
int n, m, S, T;
int deg[N];

struct node{
    int to, next;
    double c;
}e[M];
int h[N], tot;

void add(int u, int v, double c1, double c2){
    e[tot].to=v, e[tot].c=c1, e[tot].next=h[u], h[u]=tot++; 
    e[tot].to=u, e[tot].c=c2, e[tot].next=h[v], h[v]=tot++; 
}

void build(double g){
    memset(h, -1, sizeof h), tot=0;
    for(int i=1; i<=m; i++) add(edges[i].u, edges[i].v, 1, 1);
    for(int i=1; i<=n; i++) add(S, i, m, 0), add(i, T, m+2*g-deg[i], 0);
}

int q[N], d[N], cur[N];
bool bfs(){
    memset(d, -1, sizeof d);
    int tt=-1, hh=0;
    q[++tt]=S, d[S]=0, cur[S]=h[S];

    while(tt>=hh){
        int hd=q[hh++];
        for(int i=h[hd]; ~i; i=e[i].next){
            int go=e[i].to;
            if(d[go]==-1 && e[i].c>eps){
                d[go]=d[hd]+1;
                cur[go]=h[go];
                if(go==T) return true;
                q[++tt]=go;
            }
        }
    }
    return false;
}

double find(int u, double limit){
    if(u==T) return limit;
    double flow=0;
    for(int i=cur[u]; ~i && limit>flow; i=e[i].next){
        int go=e[i].to;
        cur[u]=i;
        if(d[go]==d[u]+1 && e[i].c>eps){
            double t=find(go, min(limit-flow, e[i].c));
            if(t<eps) d[go]=-1;
            e[i].c-=t, e[i^1].c+=t, flow+=t;
        }
    }
    return flow;
}

double dinic(double g){
    build(g);
    double res=0, flow;
    while(bfs()) while(flow=find(S, INF)) res+=flow;
    return res;
} 

int res=0;
bool vis[N];
void dfs(int u){
    vis[u]=true;
    if(u!=S) res++;
    for(int i=h[u]; ~i; i=e[i].next){
        int go=e[i].to;
        if(e[i].c>0 && !vis[go]) dfs(go);
    }
}

int main(){
    cin>>n>>m;

    S=0, T=n+1;
    for(int i=1; i<=m; i++){
        int u, v; cin>>u>>v;
        edges[i]={u, v};
        deg[u]++, deg[v]++;
    }

    double l=0, r=m;
    while(l+eps<r){
        double mid=(l+r)/2;
        if(m*n-dinic(mid)>eps) l=mid;
        else r=mid;
    }

    dinic(l);
    dfs(S);

    if(!res){
        puts("1\n1");
        return 0;
    }

    cout<<res<<endl;
    for(int i=1; i<=n; i++)
        if(vis[i]) cout<<i<<endl;

    return 0;
}