POJ 3155 最大密度子图

链接：

http://poj.org/problem?id=3155

题解：

要最大化下式：

可以用二分求解以下分数规划问题：

也就是最大化：

设子图为G'=(V', E')。如果边(u,v)∈E'，那么必有u,v属于V'。把点权设为负值的话，问题可以转换为求最大权闭合图(POJ 2987 Firing)。

又因为点固定时，边越多越好。所以转换思路，不是边E'决定点集V'，而应该反过来。当选定点集V'后，V'内部两两之间能构成的边就是最佳的E'。那么由V'发出的，不在E'内部的那些边构成了一个割集，当此割集最小时，E'最大。问题归结于求解最小割。

上图g为答案的猜测值，d为度数（出度+入度），c为割集（V'和V'补集构成的边）。

此时构图方式：

即是将原图 G(V , E) 转化为网络 N = (VN , EN ) 的过程：在原图点集V 的基础上增加源 s 和汇 t ；将每条原无向边 (u,v) 替换为两条容量为1的有向边u, v  和v,u  ；增加连接源 s 到原图每个点 v 的有向边s, v  ，容量为U ；增加连接原图每个点 v 到汇 t 的有向边v, t   ，容量为 (U + 2g ? dv ) 。

代码：

#include <map>
#include <set>
#include <cmath>
#include <queue>
#include <stack>
#include <cstdio>
#include <string>
#include <vector>
#include <cstdlib>
#include <cstring>
#include <sstream>
#include <iostream>
#include <algorithm>
#include <functional>
using namespace std;
#define rep(i,a,n) for (int i=a;i<n;i++)
#define per(i,a,n) for (int i=n-1;i>=a;i--)
#define pb push_back
#define mp make_pair
#define all(x) (x).begin(),(x).end()
#define SZ(x) ((int)(x).size())
typedef vector<int> VI;
typedef long long ll;
typedef pair<int, int> PII;
const ll mod = 1e9 + 7;
const int inf = 0x3f3f3f3f;
const double eps = 1e-10;
const double pi = acos(-1.0);
// head

struct edge { int to; double cap; int rev; };
const int MAX_V = 110;
vector<edge> G[MAX_V];
int level[MAX_V];
int iter[MAX_V];

void add_edge(int from, int to, double cap) {
    G[from].pb(edge{ to, cap, (int)G[to].size() });
    G[to].pb(edge{ from, 0, (int)G[from].size() - 1 });
}

void bfs(int s) {
    memset(level, -1, sizeof(level));
    queue<int> que;
    level[s] = 0;
    que.push(s);
    while (!que.empty()) {
        int v = que.front(); que.pop();
        rep(i, 0, G[v].size()) {
            edge &e = G[v][i];
            if (e.cap > 0 && level[e.to] < 0) {
                level[e.to] = level[v] + 1;
                que.push(e.to);
            }
        }
    }
}

double dfs(int v, int t, double f) {
    if (v == t) return f;
    for (int &i = iter[v]; i < G[v].size(); i++) {
        edge &e = G[v][i];
        if (e.cap > 0 && level[v] < level[e.to]) {
            double d = dfs(e.to, t, min(f, e.cap));
            if (d > 0) {
                e.cap -= d;
                G[e.to][e.rev].cap += d;
                return d;
            }
        }
    }
    return 0;
}

double max_flow(int s, int t) {
    double flow = 0;
    while (1) {
        bfs(s);
        if (level[t] < 0) return flow;
        memset(iter, 0, sizeof(iter));
        double f;
        while ((f = dfs(s, t, inf)) > 0) flow += f;
    }
}

struct Rel { int x, y; }p[1010];
int degree[MAX_V];
int n, m, s, t;

void construct_graph(double mid) {
    rep(i, 0, MAX_V) G[i].clear();
    rep(i, 1, n + 1) add_edge(s, i, m), add_edge(i, t, m + 2 * mid - degree[i]);
    rep(i, 0, m) add_edge(p[i].x, p[i].y, 1.0), add_edge(p[i].y, p[i].x, 1.0);
}

int vis[MAX_V], sum;
void dfs_travel(int v) {
    sum++;
    vis[v] = 1;
    rep(i, 0, G[v].size()){
        edge &e = G[v][i];
        if (e.cap > eps && !vis[e.to]) dfs_travel(e.to);
    }
}

int main() {
    cin >> n >> m;
    if (m == 0) return puts("1
1"), 0;
    s = 0, t = n + 1; 
    rep(i, 0, m) {
        scanf("%d%d", &p[i].x, &p[i].y);
        degree[p[i].x]++;
        degree[p[i].y]++;
    }
    double lb = 0, ub = m, precision = 1.0 / n / n;
    while (ub - lb >= precision) {
        double mid = (lb + ub) / 2;
        construct_graph(mid);
        double hg = (n*m - max_flow(s, t)) / 2;
        if (hg > eps) lb = mid;
        else ub = mid;
    }
    construct_graph(lb);
    max_flow(s, t);
    dfs_travel(0);
    cout << sum - 1 << endl;
    rep(i, 1, n + 1) if (vis[i]) cout << i << endl;
    return 0;
}