UVa 7037 ACM/ICPC 2014 Xian（网络流+最大密度子图）

7037 - The Problem Needs 3D Arrays

Time limit: 5.000 seconds

A permutation is a sequence of integers p1, p2, ..., pn, consisting of n distinct positive integers and each of them does not exceed n. Assume that r(S) of sequence S denotes the number of inversions in sequence S (if i < j and Si > Sj, then the pair of (i,j) is called an inversion of S), l(S) of sequence S denotes the length of sequence S. Given a permutation P of length n, it’s your task to find a subsequence S of P with maximum r(S) l(S) . A subsequence of P is a sequence (pi1,pi2,...,pit) which satisfies that 0 < i1 < i2 < ... < it ≤ n.
Input
The first line of the input gives the number of test cases, T. T test cases follow. For each test case, the first line contains an integer n (1 ≤ n ≤ 100), the length of the permutationP . The second line contains n integers p1, p2, ..., pn, which represents the permutation P.
Output
For each test case, output one line containing ‘Case #x: y’, where x is the test case number (starting from 1) and y is the maximum r(S) l(S) . Your answer will be considered correct if it is within an absolute error of 10−6 of the correct answer.
Sample Input
1

5

3 4 2 5 1
Sample Output
Case #1: 1.250000000000

        大致题意：给你一个1~n的排列，然后让你在这个排列中挑选任意个数字，使得这些数字组成的子串的逆序对数量除以总长最大。

        首先，我们先来回顾一下逆序对怎么计算。正常来说，O（N^2），但是大大可以用一些数据结构来优化。无非就是，加入一个数字的时候统计在此之前比该数字大的数字有多少个，然后再把这个数字加入。这种单点修改，区间查询，用树状数组再好不过了。

        但是对于这题来说，由于数字可以是任意取，即不一定是取连续的一段，所以在子串的选取上就不太好处理。如果说本题是要求子串必须是原排列的连续一段，那么直接枚举长度和其实点，用上树状数组求逆序对，O(N^2logN)即可解决问题。但是这题不是这样的，后来又尝试着用搜索加剪枝去构造子串，但是无奈还是TLE。

        后来看了看其他人的解法，网络流……的确，对于逆序对(ai,aj)，我们连边i->j，构成一个图。那么，求子串逆序对与长度之比比最大，相当于在这个图中求一个子图，使得图的密度最大，即边数除以点数最大。于是，问题就转换成了求最大密度子图的经典问题。

        说到最大密度子图，还是不得不提胡伯涛大神的那篇奇文，这个我已经引用第三次了，还是贴上来：胡伯涛《最小割模型在信息学竞赛中的应用》○|￣|_……在这里面，大神提供了两种方法。但是不管怎么样，这题首先要建立在01分数规划问题的基础上，这个最近的是在多校赛中出现过。按照通用方法，求r(S)/l(S)，我设最大值为g，那么有r(s)=l(s)*g，移项有r(s)-g*l(s)=0。然后二分g，判定结果与0的大小关系即可。但是知道了这个，我们还是需要知道怎么判定结果。这里，我们就说说两种不同的方式。

        第一种方式，我们把边也抽象为点，然后赋值上点权1，代表边数。然后原本的点保留，并赋值上点权g，对应判定式中的g*r(s)。然后，对于存在于原图中的第i条边<u,v>，我们连边<i,u>和<i,v>。代表先决条件，要取边i，那么一定要对应取点u、v。聪明的读者或许已经发现，我们要求r(s)-g*l(s)最大，可以相当于是找一个最大权闭合子图，有正负点权，有先决关系，正好符合最大权闭合子图的模型，用最小割计算即可。

        但是，我们发现，如果边比较多的时候，相应点也变多，貌似复杂度不太优，于是有了第二种方式。对于一个选定的点的集合，显然，我们选边一定要选所有的，以这个集合中点为端点的点。那么如何快速计算呢？我们考虑度的关系，如果把集合内所有的点关联的边分类，我们可以分成两类，一个是内向边，一个是外向边，图中红色为外向边，

                                                                                   

黑色加粗的为内向边。我们发现，内向边的条数可以通过度和外向边的关系求出，而外向边又可以看成集合与集合外点之间的最小割。通过以下一系列的转换，我们同样可以把这个问题转换为最小割问题。具体来说就是，首先对于所有点，源点向它连一条边，容量为U，它连向汇点，容量为U-d[i]+2*g，然后对于边<u,v>我们连u->v和v->u容量为1。其中U表示一个足够大的容量（取边的数目即可），d[i]表示i点的度。具体的证明过程见大神的文章吧，符号不好打，然后也不是很好解释。

        第二种方法，更加的复杂难懂，但是付出总是有回报，第二种方法效率非常高，一般来说比第一种方法快一倍，如果是稠密图，那么甚至可以快更多。具体见代码：

#include<bits/stdc++.h>
#define INF 0x3f3f3f3f
#define DB double
#define eps 1e-7
#define M 100010
#define N 110
using namespace std;

int n,m,a[N],d[N];
vector<int> g[N];
DB ans;

namespace ISAP
{
    int H[N],d[N],cur[N],pre[N],gap[N],Q[M];
    struct Edge{int u,v,n;DB c;} E[M];
    int nv,head,tail,cntE;DB flow,f;

    void init(){cntE=0; memset(H,-1,sizeof(H));}

    void addedge(int u,int v,DB c)
    {
        E[cntE]=Edge{u,v,H[u],c}; H[u]=cntE++;
        E[cntE]=Edge{v,u,H[v],0}; H[v]=cntE++;
    }

    void revbfs(int s,int t)
    {
        head=tail=0 ;
        memset(d,-1,sizeof(d));
        memset(gap,0,sizeof(gap));
        Q[tail++]=t;d[t]=0;gap[d[t]]=1;
        while (head!=tail)
        {
            int u=Q[head++];
            for (int i=H[u];~i;i=E[i].n)
            {
                int v=E[i].v; if (~d[v]) continue;
                d[v]=d[u]+1; gap[d[v]]++; Q[tail++]=v;
            }
        }
    }

    DB isap(int s,int t)
    {
        memcpy(cur,H,sizeof(cur)); nv=t;
        flow=0; revbfs(s,t); int u=pre[s]=s,i;
        while (d[s]<nv)
        {
            if (u==t)
            {
                f=INF;
                for (i=s;i!=t;i=E[cur[i]].v)
                    if (f-E[cur[i]].c>=eps) f=E[cur[i]].c,u=i;
                flow += f;
                for (i=s;i!=t;i=E[cur[i]].v)
                    E[cur[i]].c-=f,E[cur[i]^1].c+=f;
            }
            for (i=cur[u];~i;i=E[i].n)
                if (E[i].c>=eps&&d[u]==d[E[i].v]+1) break ;
            if (~i) cur[u]=i,pre[E[i].v]=u,u=E[i].v;
            else
            {
                if (0==--gap[d[u]]) break ;
                int minv=nv,v;
                for (int i=H[u];~i;i=E[i].n)
                {
                    v=E[i].v;
                    if (E[i].c>=eps&&minv>d[v]) minv=d[v],cur[u]=i;
                }
                d[u]=minv+1; gap[d[u]]++; u=pre[u];
            }
        }
        return flow;
    }
}

bool check(DB x)
{
    ISAP::init();
    int t=n+1,s=0;
    for(int i=1;i<=n;i++)
    {
        for(int j=0;j<g[i].size();j++)
        {
            ISAP::addedge(i,g[i][j],1);
            ISAP::addedge(g[i][j],i,1);
        }
        ISAP::addedge(s,i,m);
        ISAP::addedge(i,t,m+2*x-d[i]);
    }
    return n*m-ISAP::isap(s,t)>=eps;
}

int main()
{
    int T_T,T;
    cin>>T_T;T=T_T;
    while(T_T--)
    {
        ans=0;
        m=0; scanf("%d",&n);
        memset(g,0,sizeof(g));
        memset(d,0,sizeof(d));
        for(int i=1;i<=n;i++)
            scanf("%d",&a[i]);
        for(int i=1;i<=n;i++)
            for(int j=i+1;j<=n;j++)
                if (a[i]>a[j])
                {
                    ++m;d[i]++;d[j]++;
                    g[i].push_back(j);
                }
        DB l=0,r=(n+1)/2.0,mid;
        while(r-l>eps)
        {
            mid=(r+l)/2.0;
            if (check(mid)) ans=mid,l=mid;
                       else r=mid;
        }
        printf("Case #%d: %.6f\n",T-T_T,ans);
    }
    return 0;
}