HDU 6052 To my boyfriend 分情况处理 DP 容斥

题意：给定一个n行m列的矩阵，每个格子上有一个颜色，询问任意选择一个子矩阵，其颜色种类数的期望。n与m分范围为100。

解法：

我们定义矩阵的价值为该矩阵的颜色种类数。显然求期望的话，可以考虑求出所有子矩阵的价值之和，除以不同子矩阵的个数。

不同子矩阵的个数拿隔板法推一推可知是 (n+1)*n*(m+1)*m/4 的。

接着考虑求出所有子矩阵价值之和。

可以分别对每一种颜色进行处理。对于任一颜色，若想求出该颜色对所有子矩阵价值的贡献，可以将图看成是仅有黑白格点的图，黑点代表有颜色，白点代表无颜色，则有：贡献=总矩阵个数-全部都为白点的矩阵的个数。

求黑白格子全部都为白点的矩阵个数，涉及到一个经典题 hihocoder 1476。该题比较主流的做法有两种，一种是DP，一种是容斥：

1. DP：对于左上角点为(1, 1)，右下角点为(n, m)的矩阵，定义该矩阵内仅有白点的子矩阵其种类数为dp[i][j]，则dp[i][j]=dp[i-][j]+dp[i][j-1]-dp[i-1][j-1]+sum[以点(i, j)为右下角的全白矩阵个数]，以点(i, j)为右下角的全白矩阵个数可以用单调栈维护该矩阵的左上角可以放的位置个数。这个算法在这题中最坏的复杂度是O(nnmm)的，即当有一半的格子是黑色的时候，不过对拍了8组矩阵里只有两种颜色的数据之后，速度很快，说明常数很小。

2.容斥：全部为白点的矩阵个数=所有矩阵个数-含有一个黑点的矩阵个数+含有两个黑点的矩阵个数-含有三个黑点的矩阵个数...

含有k个黑点的矩阵个数：求出这k个黑点的上下左右边界，即能包住这k个黑点的最小矩阵，设其边界为maxx,maxy,minx,miny，含义见变量名。则总方案输为minx*miny*(n-maxx+1)*(m-maxy+1)，下标从1开始计数。复杂度是O(K*2^K)，K为黑点的个数。

但是这一题出现了一个问题：颜色种类数不定，且每种颜色的个数不定，这样子用上述两种方法都有可能返回超时(TLE)。

例如：

a) 当颜色种类数很多的时候，如每个格子颜色都不一样，这样子最多能到10000种，只用DP的做法，就会达到O(nnmm)，（这与上述那个O(nnmm)不是同一种情况），对于多组数据以及HDU的OJ环境，这样会超时；

b) 当颜色个数很多的时候，容斥方法的复杂度就指数爆炸了，也会超时。

然而，颜色的种类数和每种颜色的个数不可能同时达到最大值，故对每一种颜色，当它的点数比较少的时候，使用容斥；当它的点数比较多的时候，使用DP。

下面是AC代码：

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef pair<ll,ll> pll;
const int maxc=10005;
const int maxn=105;
int T;
bool used[maxc];
ll n,m,type;
vector<pll> col[maxc];
inline ll count_by_dp(vector<pll> &G) {
    ll dp[105][105]={0};
    for (int i=0;i<(int)G.size();++i)
        dp[G[i].first][G[i].second]=1;
    stack<pll> stk;
    int tot=0;
    for (ll i=1;i<=n;++i)
        for (ll j=1;j<=m;++j) {
            int flag=dp[i][j];
            dp[i][j]=dp[i-1][j]-dp[i-1][j-1]+dp[i][j-1];
            if (flag) {
                ++tot;
                continue;
            }
            for (int k=0;k<tot;++k) {
                if (G[k].first>i||G[k].second>j)
                    continue;
                while (!stk.empty()&&G[k].second>=stk.top().second)
                    stk.pop();
                stk.push(G[k]);
            }
            ll sum=i*j,prey=0;
            while (!stk.empty()) {
                ll x=stk.top().first;
                ll y=stk.top().second;
                stk.pop();
                sum-=x*(y-prey);
                prey=y;
            }
            dp[i][j]+=sum;
        }
    return dp[n][m];
}
inline ll count_by_inclusion_exclusion(vector<pll> &G) {
    int K=G.size(),sz=1<<K;
    ll ret=(n+1)*n*(m+1)*m/4LL;
    for (int i=1;i<sz;++i) {
        int coe=1;
        ll minx=n,miny=m;
        ll maxx=0,maxy=0;
        for (int j=0;j<K;++j)
            if (i&(1<<j)) {
                minx=min(minx,G[j].first);
                miny=min(miny,G[j].second);
                maxx=max(maxx,G[j].first);
                maxy=max(maxy,G[j].second);
                coe=-coe;
            }
        ret+=coe*minx*miny*(n-maxx+1)*(m-maxy+1);
    }
    return ret;
}
int main()
{
    scanf("%d",&T);
    while (T--) {
        memset(used,0,sizeof used);
        type=0;
        scanf("%I64d%I64d",&n,&m);
        for (int i=0;i<n*m;++i)
            col[i].clear();
        for (int i=1;i<=n;++i)
            for (int j=1;j<=m;++j) {
                int a;
                scanf("%d",&a);
                col[a].push_back(pll(i,j));
                if (!used[a]) {
                    ++type;
                    used[a]=1;
                }
            }
        ll res=type*n*(n+1)*m*(m+1)/4LL;
        ll down=n*(n+1)*m*(m+1)/4LL;
        for (int i=0;i<n*m;++i) {
            sort(col[i].begin(),col[i].end());
            int sz=col[i].size();
            if (sz>13)
                res-=count_by_dp(col[i]);
            else if (sz)
                res-=count_by_inclusion_exclusion(col[i]);
        }
        printf("%.9f\n",1.0*res/down);
    }
    return 0;
}