[TC-HouseProtection]House Protection

题目大意：
　　一个平面直角坐标系中有给定的$n(n\le50)$个红点和$m(m\le50)$个蓝点，每个点可以选择画一个半径为$r$（所有的$r$相同）的圆或不画。圆的半径上限为$R(R\le1000)$。且不同颜色的点所画成的圆不能相交，问所有圆面积的和最大是多少？

思路：
　　枚举每一对不同颜色的点，求出所有可能的$r$，显然这样的$r$有$nm$个。
　　对于不同颜色的点，若以$r$为半径画出的圆相交，则在这两个点之间连边，题目转化为最大独立集问题。有结论：最大独立集=点数-二分图最大匹配。用Dinic求最大匹配即可。
　　对$r$排序，则每次只会增加新的边，只需要在原来的基础上进行增广即可。时间复杂度$O(n^4)$。

#include<cmath>
#include<queue>
#include<vector>
#include<cstring>
#include<climits>
#include<algorithm>
class HouseProtection {
    private:
        static constexpr int N=50,M=102;
        static constexpr double pi=M_PI;
        using Point=std::pair<int,int>;
        struct Length {
            double l;
            int u,v;
            bool operator < (const Length &rhs) const {
                return l<rhs.l;
            }
        };
        struct Edge {
            int from,to,remain,next;
        };
        Point a[N],b[N];
        std::vector<Length> len;
        std::vector<Edge> e;
        int s,t,lev[M],cur[M],h[M];
        double sqr(const double &x) const {
            return x*x;
        }
        double dist(const Point &a,const Point &b) const {
            return sqrt(sqr(a.first-b.first)+sqr(a.second-b.second));
        }
        void add_edge(const int &u,const int &v,const int &w) {
            e.push_back({u,v,w,h[u]});
            h[u]=e.size()-1;
        }
        void bfs() {
            memset(lev,-1,sizeof lev);
            lev[s]=0;
            std::queue<int> q;
            q.push(s);
            while(!q.empty()) {
                const int &x=q.front();
                for(register int i=h[x];~i;i=e[i].next) {
                    const int &y=e[i].to,&r=e[i].remain;
                    if(r&&!~lev[y]) {
                        lev[y]=lev[x]+1;
                        q.push(y);
                    }
                }
                q.pop();
            }
        }
        int dfs(const int &x,const int &flow) {
            if(x==t) return flow;
            for(int &i=cur[x];~i;i=e[i].next) {
                const int &y=e[i].to;
                if(e[i].remain&&lev[y]>lev[x]) {
                    if(int d=dfs(y,std::min(flow,e[i].remain))) {
                        e[i].remain-=d;
                        e[i^1].remain+=d;
                        return d;
                    }
                }
            }
            return 0;
        }
        int dinic() {
            int maxflow=0;
            for(;;) {
                bfs();
                if(!~lev[t]) break;
                memcpy(cur,h,sizeof h);
                while(int flow=dfs(s,INT_MAX)) {
                    maxflow+=flow;
                }
            }
            return maxflow;
        }
    public:
        double safetyFactor(const std::vector<int> &possibleXForBlue,const std::vector<int> &possibleYForBlue,const std::vector<int> &possibleXForRed,const std::vector<int> &possibleYForRed,const int &R) {
            memset(h,-1,sizeof h);
            const int n=possibleXForBlue.size(),m=possibleXForRed.size();
            for(register int i=0;i<n;i++) {
                a[i]={possibleXForBlue[i],possibleYForBlue[i]};
            }
            for(register int i=0;i<m;i++) {
                b[i]={possibleXForRed[i],possibleYForRed[i]};
            }
            for(register int i=0;i<n;i++) {
                for(register int j=0;j<m;j++) {
                    const double dis=dist(a[i],b[j])/2;
                    if(dis<R) len.push_back({dis,i,j});
                }
            }
            len.push_back({(double)R,-1,-1});
            std::sort(len.begin(),len.end());
            s=n+m,t=n+m+1;
            for(int i=0;i<n;i++) {
                add_edge(s,i,1);
                add_edge(i,s,0);
            }
            for(int i=0;i<m;i++) {
                add_edge(n+i,t,1);
                add_edge(t,n+i,0);
            }
            int match=0;
            double ans=0;
            for(register auto &e:len) {
                match+=dinic();
                ans=std::max(ans,pi*sqr(e.l)*(n+m-match));
                if(e.l<R) {
                    add_edge(e.u,n+e.v,1);
                    add_edge(n+e.v,e.u,0);
                }
            }
            return ans;
        }
};