【GDOI2017模拟11.2】抗拒黄泉

原题：topcoder SRM583 的 Hard

题目大意

有n*m的棋盘，每个格子的值为0或1，每次我可以等概率地选择一个1，并给这个位置打上标记（一个格子可以重复标记）。求：当每一行、每一列都有至少一个标记时，标记次数的期望。

n,m≤20 n*m≤200

分析

设P(i)表示进行i次操作后没有到达最终状态的概率，那么答案Ans=$\sum_{i=0}^\infty P(i)$

现在考虑如何计算它。假设经过了i次，我有一个行列状态S，这个状态表示的行与列都没有被选取，p(S)表示一次选中它的概率。由于我i次都没有选中它，概率就为$(1-p(S))^i$。但是这样算不能确定其它行列是否被选，就要用容斥。那么

$$P(i)=\sum_{|S|为奇数}(1-p(S))^i- \sum_{|S|为偶数}(1-p(S))^i$$
那么 $$Ans=\sum_{i=0}^\infty (\sum_{|S|为奇数}(1-p(S))^i- \sum_{|S|为偶数}(1-p(S))^i)$$
括号里面的东西看成两个等比数列
化简，得 $$Ans=\sum_{|S|为奇数}\frac{1}{p(S)}- \sum_{|S|为偶数}\frac{1}{p(S)}$$
设sum(S)表示组成行列状态S的1的个数（sum(S)可以为很多值），tot为1的总数，那么对于一种可能，p(S)=sum(S)/tot，那么 $$Ans=\sum_{|S|为奇数}\frac{tot}{sum(S)}- \sum_{|S|为偶数}\frac{tot}{sum(S)}$$
如果确定了S，可以枚举sum(S)的值，求有多少个可能。由于行+列较大，可以考虑枚举行的状态（这样就锁定了一些行的状态），然后做一个dp：设f[i][j][k]表示当前确定了前i列的状态，|S|的奇偶性为j，元素个数为k，的方案数。这只是一个背包。

时间复杂度$O(2^n*m^2n)$

观察到n*m≤200，把n和m中较小的作为行，这样跑得快一些。

#include <cstdio>
#include <cstring>
#include <algorithm>

using namespace std;

const int maxn=205;

typedef long long LL;

typedef double db;

int n,m,N,M,t,map[25][25],c[25][25],tot,f[25][2][maxn];

bool bz[25];

double ans;

void dp()
{
    int t=0,sum=0;
    for (int i=1;i<=n;i++) if (bz[i])
    {
        t^=1;
        for (int j=1;j<=m;j++) sum+=map[i][j];
    }
    memset(f,0,sizeof(f));
    f[0][t][sum]=1;
    for (int i=1;i<=m;i++)
    {
        int cnt=0;
        for (int j=1;j<=n;j++) if (!bz[j] && map[j][i]>0) cnt++;
        for (int j=0;j<2;j++)
        {
            for (int k=0;k<=tot-cnt;k++) if (f[i-1][j][k])
            {
                f[i][j][k]+=f[i-1][j][k];
                f[i][j^1][k+cnt]+=f[i-1][j][k];
            }
        }
    }
    for (int i=1;i<=tot;i++)
    {
        ans+=(db)f[m][1][i]*tot/i-(db)f[m][0][i]*tot/i;
    }
}

void dfs(int x)
{
    if (x>n)
    {
        dp(); return;
    }
    bz[x]=0; dfs(x+1);
    bz[x]=1; dfs(x+1);
}

int main()
{
    freopen("refuse.in","r",stdin); //freopen("refuse.out","w",stdout);
    scanf("%d%d",&n,&m);
    for (int i=1;i<=n;i++) for (int j=1;j<=m;j++)
    {
        scanf("%d",&c[i][j]);
        tot+=c[i][j];
    }
    if (n>m)
    {
        for (int i=1;i<=m;i++)
        {
            for (int j=1;j<=n;j++) map[i][j]=c[j][i];
        }
        n^=m^=n^=m;
    }else memcpy(map,c,sizeof(map));
    dfs(1);
    printf("%.5lf\n",ans);
    fclose(stdin); fclose(stdout);
    return 0;
}