poj 3155 01规划->最大密度密度子图->最大流

对于求最大密度子图，就是求

max{|E|/|V|}

然后就可以写为a(x)为所选的边，b(x)为所选的点。（边显然是点的导出子图）

然后a(x)>0,b(x)>0，也是显然的条件。

然后就可以改为和01规划一样，求 f(g) = max{a(x)-g*b(x)} ，这个函数是单调递减的。当f(g)=0能取得最优解（证明略）

然后这都是论文的内容。。

然后就是构图。

s->每个点，流量为巨大值，不是无限大。。

原来每个点之间，正反两条边，流量都是1

每个点->t ,流量为2g-dv+巨大值  ，其中g为f(g)的g,dv为每个点的度。

巨大值=边的总数。

然后根据论文的性质，任意两个方案，他们结果的差再1/(n^2)之间，这也就决定了二分答案找g的精度。

然后当g大的时候，f(g)为负数。但是找最大流的时候，会出现“0”流量的情况，这要进行考虑，要特判。

然后差不多就做完了。

再次证明了，板子速度还是蛮快的嘛

ac code

#include <bits/stdc++.h>
#include <ext/pb_ds/priority_queue.hpp>
#include <tr1/unordered_map>
using std::tr1::unordered_map;
//using std::setiosflags;
//using std::setprecision;
using std::sort;
using std::max;
using std::min;
using std::cout;
using std::stack;
using std::cin;
using std::endl;
using std::swap;
using std::pair;
using std::vector;
using std::set;
using std::map;
using std::make_pair;
using std::multiset;
using std::unique;
using std::queue;
using std::greater;
using std::string;
using std::priority_queue;
using std::lower_bound;//返回第一个不小于
using std::upper_bound;//返回第一个大于
using std::max_element;
using std::min_element;
using __gnu_pbds::pairing_heap_tag;
#define x first
#define y second
#define Hash unordered_map
#define clr(x,b) memset((x),(b),sizeof(x))
typedef unsigned long long uLL;
typedef long long LL;
typedef pair<int, int> pii;
typedef pair<double, double> pdd;
typedef __gnu_pbds::priority_queue<pii, greater<pii>, pairing_heap_tag> Heap;//小根堆
typedef Heap::point_iterator Hit;
const Hit null;
const double PI = acos(-1);
const LL LINF = 0x3f3f3f3f3f3f3f3fll;//4e18
//const int INF = 0x3f3f3f3f;//1e9
const double INF = 1e13;
const double eps = 1e-7;
const int MOD = 1 << 16;
#define prln(x) cout<<#x<<" = "<<x<<endl
#define pr(x) cout<<#x<<" = "<<x<<" "

//调用方法：
//init(n) 初始化 n 为节点数
//clear_flow() 清空所有边的流量
//add_edge(from, to, cap) 添加边from -> to 容量为cap
//maxflow(s, t) 返回s->t的最大流
//void mincut(vector<int>& ans) // 调用完maxflow后才可以用，ans里面存最小割
//print() 打印整张图，调试用
const int maxn = 200;

template<typename T>
struct Edge { int from, to; T cap, flow;};

template<typename T>
struct ISAP {
int n, m, s, t;
vector<Edge<T> > edges;
vector<int> G[maxn]; // 邻接表，G[i][j]表示结点i的第j条边在e数组中的序号
bool vis[maxn]; // BFS使用
int d[maxn]; // 从起点到i的距离
int cur[maxn]; // 当前弧指针
int p[maxn]; // 可增广路上的上一条弧
int num[maxn]; // 距离标号计数
void add_edge(int from, int to, T cap) {
//pr(from),pr(to),prln(cap);
edges.push_back((Edge<T>){from, to, cap, 0});
edges.push_back((Edge<T>){to, from, 0, 0});
m = edges.size();
G[from].push_back(m-2);
G[to].push_back(m-1);
}
bool bfs() {
memset(vis, 0, sizeof(vis));
queue<int> q;
q.push(t);
vis[t] = 1;
d[t] = 0;
while(!q.empty()) {
int x = q.front(); q.pop();
for(int i = 0; i < G[x].size(); i++) {
Edge<T>& e = edges[G[x][i]^1];
if(!vis[e.from] && e.cap > e.flow) {
vis[e.from] = 1;
d[e.from] = d[x] + 1;
q.push(e.from);
}
}
}
return vis[s];
}

vector<int>output;
queue<int>q;
void solve()
{
clr(vis,0);
vis[s]=1;
output.clear();
q.push(s);
while (!q.empty())
{
int x = q.front(); q.pop();
//prln(x);
for (int i = 0; i < G[x].size(); i ++)
{
Edge<T>& e = edges[G[x][i]];
if (!vis[e.to] && e.cap > e.flow)
{
vis[e.to] = 1;
q.push(e.to);
output.push_back(e.to);
}
}
}
printf("%d\n",output.size());
sort(output.begin(), output.end());
for (auto x : output)
{
printf("%d\n",x);
}
}

void init(int n) {
this->n = n;
for(int i = 0; i < n; i++) G[i].clear();
edges.clear();
}
void clear_flow() {
for(int i = 0; i < edges.size(); i++) edges[i].flow = 0;
}
double augment() {
int x = t;
T a = INF;//初始增广流量
while(x != s) {
Edge<T>& e = edges[p[x]];
a = min(a, e.cap-e.flow);
x = edges[p[x]].from;
}
x = t;
while(x != s) {
edges[p[x]].flow += a;
edges[p[x]^1].flow -= a;
x = edges[p[x]].from;
}
return a;
}
T maxflow(int s, int t) {//到的最大流大于need就停止，如果没有限制，删去含有need的地方
this->s = s; this->t = t;
T flow = 0;
bfs();
memset(num, 0, sizeof(num));
for(int i = 0; i < n; i++) num[d[i]]++;
int x = s;
memset(cur, 0, sizeof(cur));
while(d[s] < n) {
if(x == t) {
flow += augment();
x = s;
}
int ok = 0;
for(int i = cur[x]; i < G[x].size(); i++) {
Edge<T>& e = edges[G[x][i]];
if(e.cap > e.flow && d[x] == d[e.to] + 1) { // Advance
ok = 1;
p[e.to] = G[x][i];
cur[x] = i; // 注意
x = e.to;
break;
}
}
if(!ok) { // Retreat
int m = n-1; // 初值注意
for(int i = 0; i < G[x].size(); i++) {
Edge<T>& e = edges[G[x][i]];
if(e.cap > e.flow) m = min(m, d[e.to]);
}
if(--num[d[x]] == 0) break;//gap优化
num[d[x] = m+1]++;
cur[x] = 0; // 注意
if(x != s) x = edges[p[x]].from;
}
}

//判断是否是一个解，如果是所有点都没有选中，就是非法的
for (int i = 0; i < G[s].size(); i ++ )
{
Edge<T>& e = edges[G[s][i]];
if (e.cap > e.flow)//如果已经有边流量可以出，就是合法解。非法解的话，直接干掉
{
return flow;
}
}
return INF;
}

void mincut(vector<int>& ans) { // 调用完maxflow后才可以用，ans里面存最小割
bfs();
for(int i = 0; i < edges.size(); i++) {
Edge<T>& e = edges[i];
if(!vis[e.from] && vis[e.to] && e.cap > 0) ans.push_back(i);
}
}

void print() {
printf("Graph:\n");
for(int i = 0; i < edges.size(); i++)
printf("%d->%d, %.5lf, %.5lf\n", edges[i].from, edges[i].to , edges[i].cap, edges[i].flow);
}
} ;

ISAP<double>isap;

int n, m;
struct bian
{
int from, to;
}eee[2000];

int degree[maxn];

//最大密度子图，密度不超过m/1，不小于1/n。
//任意两个密度之间，答案误差不小于1/n^2。 当精度到这个误差内的时候，就可以结束了
void init()
{
isap.init(n + 5);
clr(degree, 0);
for (int i = 0; i < m; ++ i)
{
scanf("%d%d", &eee[i].from, &eee[i].to);
++ degree[eee[i].from];
++ degree[eee[i].to];
}
}

int U;

double check(double lamuda)
{
//求min{lamuda * b(x) - a(x)},其中b(x)为点，a(x)为边
//求max{a(x) - lamuda * b(x)},两个显然是等价的
//根据论文构图，最终结果
//最大流为最小割，最小割是U*n-2(|E'|-g|V'|)的最小值
//此时，显然是|E'|-g|V'|的最大值。
//也就是用(U * n-maxflow)/2，为|E'|-g|V'|的最大值，正好是规划的要求。查看是否为0即可.
double U = m;//即可保证2*lamuda-dv为为正数。这个式子见论文。
isap.init(n + 5);
for (int i = 1; i <= n; ++ i)
{
//这些边的编号都为0
isap.add_edge(0,i,U);
isap.add_edge(i, n+1, U + 2 * lamuda - degree[i]);
}
for (int i = 0; i < m; ++ i)
{
isap.add_edge(eee[i].from, eee[i].to, 1);
isap.add_edge(eee[i].to, eee[i].from, 1);
}
double ret = isap.maxflow(0, n + 1);
//prln(ret);
//printf("%.6lf %.6lf\n",ret, U*n);
return (U * n - ret);
}

void doit()
{
double R = 1.0*m/1, L = 1.0*1/n, mid;

while (R-L >= + 1.0/(n*n))//精度在一定范围内即可
{
//printf("%.6lf %.6lf %.6lf\n", L, R, check((L+R)/2));
mid = L + (R-L)/2;
if (check(mid) < 0) R = mid;
else L = mid;
}
//R和mid都可能是非法解，L一定是合法解
check(L);
}

int main()
{
// freopen("hard.in","r",stdin);
// freopen("hard.out","w",stdout);
while (~scanf("%d%d", &n, &m))
{
if (m==0)
{
cout<<1<<endl<<1<<endl;
continue;
}
init();
doit();
isap.solve();
}
return 0;
}