【网络流】最大密度子图

参考：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_{uin V'}d_u - cnt[V',overline{V'}]}{2} ]

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

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

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

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

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

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

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

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

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

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

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

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

进一步化为：

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

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

[c[s,t]=|V|U + sum_{uin 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
1");
        return 0;
    }

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

    return 0;
}