[网络流24题] 04 魔术球问题 (有向无环图最小路径覆盖, 最大流)

题目大意：

假设有n根这柱子， 现在按照下述规则在这n根柱子中一次放入编号为1， 2， 3， 4 .......的球。

①：每次只能在某根柱子的最上面放球。

②：在同一根柱子中， 任意两个相邻球的编号之和为完全平方数。

现在给出n根柱子， 计算最多能放多少个球。

思路分析：

枚举答案转化为判定性问题，然后最小路径覆盖，可以转化成二分图最大匹配，从而用最大流解决。

枚举答案a，1，2，3， .... ,a，在图中建立结点， 设立一个源点一个汇点；

<建模方法>：

①：每次新增一个球cnt， 属于X集(cnt)让其与源点相连， 对应的Y集元素(cnt+NN)与汇点相连， 容量都为1；

②：让X集的元素，不包括刚刚新增的球，每一个元素满足与Y集刚刚新增的cnt+NN元素的和为完全平方数， 则从X集该元素向Y集的新增的球连一条边， 容量为1；

③：不断的跑最大流， 求最小路径覆盖。如果刚好满足最小路径覆盖数等于n，a-1就是最优解；

<建模分析>：

由于是顺序放球，每根柱子上的球满足这样的特征，即下面的球编号小于上面球的编号。抽象成图论，把每个球看作一个顶点，就是编号较小的顶点向编号较大的顶点连接边，条件是两个球可以相邻，即编号之和为完全平方数。每根柱子看做一条路径，n根柱子要覆盖掉所有点，一个解就是一个路径覆盖。最小路径覆盖数随球的数量递增不递减，满足单调性，所以可以枚举答案（或二分答案），对于特定的答案求出最小路径覆盖数，一个可行解就是最小路径覆盖数等于n的答案，求出最大的可行解就是最优解。本问题更适合枚举答案而不是二分答案，因为如果顺序枚举答案，每次只需要在残量网络上增加新的节点和边，再增广一次即可。如果二分答案，就需要每次重新建图，大大增加了时间复杂度。

代码实现：（注释部分为输出路径）

#include<cstdio>
#include<cstring>
#include<queue>
#include<iostream>
#include<algorithm>
#include<cmath>
using namespace std;
const int N=4001, M=200000, NN=2000, INF=0x3f3f3f3f;
int n, m, s, t, top, cnt;
int head[N], dis[N], gap[N], pre[N], cur[N], path[N];
bool vis[N], is_sqr[N];

struct Edge{
    int to, next, flow;
    Edge(int _to = 0, int _next = 0, int _flow = 0):to(_to), next(_next), flow(_flow){}
}edge[M];

void Addedge(int from, int to, int flow){
    edge[top] = Edge(to, head[from], flow);
    head[from] = top++;
    edge[top] = Edge(from, head[to], 0);
    head[to] = top++;
}

void Bfs(){
    queue<int> q;
    memset(gap, 0, sizeof(gap));
    memset(dis, -1, sizeof(dis));
    gap[0] = 1, dis[t] = 0;
    q.push(t);
    while(!q.empty()){
        int u = q.front(); q.pop();
        for(int i = head[u]; i + 1; i = edge[i].next){
            if(dis[edge[i].to] == -1){
                dis[edge[i].to] = dis[u] + 1;
                gap[dis[edge[i].to]] ++;
                q.push(edge[i].to);
            }
        }
    }
}

int Sap(){
    Bfs();
    memset(pre, -1, sizeof(pre));
    for(int i = 0; i <= t; ++i) cur[i] = head[i];
    int u = s, cur_flow, max_flow = 0, neck, tmp;
    while(dis[s] <= t){
        if(u == t){
            cur_flow = INF;
            for(int i = s; i != t; i = edge[cur[i]].to){
                if(cur_flow > edge[cur[i]].flow){
                    neck = i;
                    cur_flow = edge[cur[i]].flow;
                }
            }
            for(int i = s; i != t; i = edge[cur[i]].to){
                tmp = cur[i];
                edge[tmp].flow -= cur_flow;
                edge[tmp^1].flow += cur_flow;
            }
            max_flow += cur_flow;
            u = neck;
        }
        int i;
        for(i = cur[u]; i + 1; i = edge[i].next)
            if(edge[i].flow && dis[u] == dis[edge[i].to] + 1) break;
        if(i != -1){
            cur[u] = i;
            pre[edge[i].to] = u;
            u = edge[i].to;
        }else{
            if(--gap[dis[u]] == 0) break;
            cur[u] = head[u];
            int mindis = t;
            for(i = head[u]; i + 1; i = edge[i].next){
                if(edge[i].flow && mindis > dis[edge[i].to]) mindis = dis[edge[i].to];
            }
            dis[u] = mindis + 1;
            gap[dis[u]]++;
            if(u != s) u = pre[u];
        }
    }
    return max_flow;
}

//void Dfs(int u){
//    vis[u] = 1;
//    path[cnt++] = u;
//    for(int i = head[u]; i+1; i = edge[i].next){
//        if(!vis[i] && edge[i].to && !edge[i].flow) Dfs(edge[i].to - NN);
//    }
//}

int main(){
    freopen("balla.in", "r", stdin);
    freopen("balla.out", "w", stdout);
    memset(head, -1, sizeof(head));
    memset(is_sqr, 0, sizeof(is_sqr));
    memset(vis, 0, sizeof(vis));
    scanf("%d", &n);
    s = top = cnt = 0, t = N - 1;
    int res = 0, ans = 0;
    for(int i = 1; i < 60; ++i) is_sqr[i*i] = 1;
    while(1){
        cnt++;
        Addedge(s, cnt, 1), Addedge(cnt+NN, t, 1);
        for(int i = 1; i < cnt; ++i) if(is_sqr[i+cnt]) Addedge(i, cnt+NN, 1);
        res += Sap();
        if(cnt - res > n) break;
    }
    ans = cnt-1;
    printf("%d\n", ans);
//    memset(head, -1, sizeof(head));
//    top = 0;
//    for(int i = 1; i <= ans; ++i){
//        Addedge(s, i, 1), Addedge(i+NN, t, 1);
//        for(int j = 1; j < i; ++j)
//            if(is_sqr[i+j]) Addedge(j, i+NN, 1);
//    }
//    res = Sap();
//    for(int i = 1; i <= ans; ++i){
//        cnt = 0;
//        if(!vis[i]){
//            Dfs(i);
//            printf("%d", path[0]);
//            for(int j = 1; j < cnt; ++j) printf(" %d", path[j]);
//            printf("\n");
//        }
//    }
}