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 ) 。
代码：
1 #include <map> 2 #include <set> 3 #include <cmath> 4 #include <queue> 5 #include <stack> 6 #include <cstdio> 7 #include <string> 8 #include <vector> 9 #include <cstdlib> 10 #include <cstring> 11 #include <sstream> 12 #include <iostream> 13 #include <algorithm> 14 #include <functional> 15 using namespace std; 16 #define rep(i,a,n) for (int i=a;i<n;i++) 17 #define per(i,a,n) for (int i=n-1;i>=a;i--) 18 #define pb push_back 19 #define mp make_pair 20 #define all(x) (x).begin(),(x).end() 21 #define SZ(x) ((int)(x).size()) 22 typedef vector<int> VI; 23 typedef long long ll; 24 typedef pair<int, int> PII; 25 const ll mod = 1e9 + 7; 26 const int inf = 0x3f3f3f3f; 27 const double eps = 1e-10; 28 const double pi = acos(-1.0); 29 // head 30 31 struct edge { int to; double cap; int rev; }; 32 const int MAX_V = 110; 33 vector<edge> G[MAX_V]; 34 int level[MAX_V]; 35 int iter[MAX_V]; 36 37 void add_edge(int from, int to, double cap) { 38 G[from].pb(edge{ to, cap, (int)G[to].size() }); 39 G[to].pb(edge{ from, 0, (int)G[from].size() - 1 }); 40 } 41 42 void bfs(int s) { 43 memset(level, -1, sizeof(level)); 44 queue<int> que; 45 level[s] = 0; 46 que.push(s); 47 while (!que.empty()) { 48 int v = que.front(); que.pop(); 49 rep(i, 0, G[v].size()) { 50 edge &e = G[v][i]; 51 if (e.cap > 0 && level[e.to] < 0) { 52 level[e.to] = level[v] + 1; 53 que.push(e.to); 54 } 55 } 56 } 57 } 58 59 double dfs(int v, int t, double f) { 60 if (v == t) return f; 61 for (int &i = iter[v]; i < G[v].size(); i++) { 62 edge &e = G[v][i]; 63 if (e.cap > 0 && level[v] < level[e.to]) { 64 double d = dfs(e.to, t, min(f, e.cap)); 65 if (d > 0) { 66 e.cap -= d; 67 G[e.to][e.rev].cap += d; 68 return d; 69 } 70 } 71 } 72 return 0; 73 } 74 75 double max_flow(int s, int t) { 76 double flow = 0; 77 while (1) { 78 bfs(s); 79 if (level[t] < 0) return flow; 80 memset(iter, 0, sizeof(iter)); 81 double f; 82 while ((f = dfs(s, t, inf)) > 0) flow += f; 83 } 84 } 85 86 struct Rel { int x, y; }p[1010]; 87 int degree[MAX_V]; 88 int n, m, s, t; 89 90 void construct_graph(double mid) { 91 rep(i, 0, MAX_V) G[i].clear(); 92 rep(i, 1, n + 1) add_edge(s, i, m), add_edge(i, t, m + 2 * mid - degree[i]); 93 rep(i, 0, m) add_edge(p[i].x, p[i].y, 1.0), add_edge(p[i].y, p[i].x, 1.0); 94 } 95 96 int vis[MAX_V], sum; 97 void dfs_travel(int v) { 98 sum++; 99 vis[v] = 1; 100 rep(i, 0, G[v].size()){ 101 edge &e = G[v][i]; 102 if (e.cap > eps && !vis[e.to]) dfs_travel(e.to); 103 } 104 } 105 106 int main() { 107 cin >> n >> m; 108 if (m == 0) return puts("1\n1"), 0; 109 s = 0, t = n + 1; 110 rep(i, 0, m) { 111 scanf("%d%d", &p[i].x, &p[i].y); 112 degree[p[i].x]++; 113 degree[p[i].y]++; 114 } 115 double lb = 0, ub = m, precision = 1.0 / n / n; 116 while (ub - lb >= precision) { 117 double mid = (lb + ub) / 2; 118 construct_graph(mid); 119 double hg = (n*m - max_flow(s, t)) / 2; 120 if (hg > eps) lb = mid; 121 else ub = mid; 122 } 123 construct_graph(lb); 124 max_flow(s, t); 125 dfs_travel(0); 126 cout << sum - 1 << endl; 127 rep(i, 1, n + 1) if (vis[i]) cout << i << endl; 128 return 0; 129 }
