【2017集训队作业】小c的岛屿（斯特林反演）（期望）（高斯消元）

题面：

有 $n$ 座岛屿，初始时没有边。每座岛屿都有一个概率值 $p_i$ 和一个友好列表 $A_i$。

小 c 站在 $1$ 号岛屿，依次执行以下操作：

1. 设现在在岛屿 $x$，有 $p_x$ 的概率产生一条图中不存在的随机无向边，不会产生自环。
2. 如果此时所有岛屿仍未连通，他会在当前点的友好列表中，随机选择一个，走到那座岛屿上，并把不满意度 $+1$，然后重复第 $1$ 步。否则结束这个过程。

求他的期望不满意度，模 $10^9+7$。
$n\le 50$。

---

题解：

注意1步骤中产生的边与2步骤中的友好列表无关，就算两点之间没有边，只要友好列表中存在$u\to v$，$u$就能直接走到$v$。

感觉很妙妙的一道题，记录一下思路。

O(n^5)解法

设$E(u,i)$表示当前处于点$u$，已经连了$i$条边，之后能够获得的期望不满意度。

根据题目我们知道所有$n$个点$m$条边的图产生概率是相等的。

很容易地，我们可以得到方程：$$E(u,i)=p_u\cdot h_{i+1}\cdot \frac{1}{d_u}\sum _{u\to v}(E(v,i+1)+1)+(1-p_u)\frac{1}{d_u}\sum _{u\to v}(E(v,i)+1)$$

其中$d_u$表示点$u$的度数，$h_i$表示$n$个点$i-1$条边的非连通图，加一条边之后仍然不联通的概率。

很显然上面这个东西考虑分层高斯消元就是$O(n^5)$，常数优秀是能过的，现在的问题在于怎么求$h$。

显然对于$h$，有$$h_{i+1}=\frac{i+1条边非连通图个数\times (i+1)}{(C_n^2-i)\times i条边非连通图个数}$$

问题变成了求$n$个点若干条边的连通图个数。

显然有一种枚举一号点所在的连通块的DP做法，复杂度$O(n^6)$不太行。

考虑容斥，计算连通块划分后斯特林反演。

设考虑将$n$个点划分为$k$个连通块，不同块之间不允许连边，同一块内可以任意连边（哪怕连出来不同）。则一个实际连通块数量为$t$的图会在计算$k$个连通块的时候算到$S_{t,k}$次。乘上$(-1{)}^{k-1}(k-1)!$的系数直接斯特林反演得到一个连通块的情况。

不过别忘了，我们这道题有关于边的变量。

所以$i$条边，$n$个点的连通图数量实际上为：

$$\sum _G(-1{)}^{c(G)-1}(c(G)-1)!(\genfrac{}{}{0ex}{}{e(G)}{i})$$

其中$G$需要取遍所有的$n$个点的，划分成了若干连通块的，每个连通块内连成了完全图的图，$c(G)$表示$G$的连通块数量，$e(G)$表示$G$的边数。

设$f[i][j]$表示$i$个点，划分成了若干连通块，每个连通块内连成完全图后边数和为$j$，的所有图的$(-1{)}^{c(G)-1}(c(G)-1)!$之和。

感觉这个$(c(G)-1)!$并不是很好处理。

我们可以先定一号店所在的连通块大小，然后转移的时候枚举下一个连通块大小，乘上$-1$转移，这样一个$c(G)$个连通块的图刚好被算了$(c(G)-1)!$次。

这样可以$O(n^4)$计算出$f[n][i]$。然后$O(n^4)$直接算出$g[i]$，表示$n$个点$i$条边的连通图个数。

$h_i$可以直接计算，然后就是分层高斯消元了。把变量设置成$(E(u,i)+1)$来高斯消元会好写点。

O(n^4)解法

上面算法的复杂度瓶颈在于分层的高斯消元。

而分层的原因是每层的$h_i$不同。

于是换一种思路，设$E(u,i)$表示当前在点$u$，需要再得到$i$条边的期望步数。

则答案为$$\sum E(1,i)\cdot P(i条边恰好使图连通的概率)$$

而$P(i条边恰好使图连通的概率)=P(i条边的图连通概率)-P(i-1条边的图连通概率)$

算一下连通图的数量就行了。

然后我们知道有方程$$E(u,i)=1+\frac{1}{d_u}\sum _{u\to v}{p}_{u}E(v,i-1)+(1-p_u)E(v,i)$$

递推系数是固定的，高斯消元解一下递推系数就行了。

然后就可以$O(n^4)$递推了

---

代码（复杂度$O(n^5)$）：

#include<bits/stdc++.h>
#define ll long long
#define re register
#define cs const

using std::cerr;
using std::cout;

cs int mod=1e9+7;
inline int add(int a,int b){a+=b-mod;return a+((a>>31)&mod);}
inline int dec(int a,int b){a-=b;return a+((a>>31)&mod);}
inline int mul(int a,int b){ll r=(ll)a*b;return r>=mod?r%mod:r;}
inline int power(int a,int b,int res=1){
	for(;b;b>>=1,a=mul(a,a))(b&1)&&(res=mul(res,a));
	return res;
}
inline void Inc(int &a,int b){a+=b-mod;a+=(a>>31)&mod;}
inline void Dec(int &a,int b){a-=b;a+=(a>>31)&mod;}
inline void Mul(int &a,int b){a=mul(a,b);} 

cs int N=55,M=N*N>>1;

int C[M][M];
int c[N];
inline void init(){
	for(int re i=0;i<M;++i)C[i][0]=1;
	for(int re i=1;i<M;++i)
	for(int re j=1;j<=i;++j)C[i][j]=add(C[i-1][j],C[i-1][j-1]);
	for(int re i=1;i<N;++i)c[i]=i*(i-1)>>1;
}

int n,m;
std::vector<int> G[N];
int p[N],f[M][M],g[M],h[M];

int A[N][N];
inline void gauss(int mm,int *a,int *b){
	for(int re i=1;i<=n;++i)memset(A[i]+1,0,sizeof(int)*(n+1));
	for(int re i=1;i<=n;++i){
		A[i][i]=mod-1;
		Inc(A[i][n+1],1);
		int inv=power(G[i].size(),mod-2);
		for(int re v:G[i]){
			Inc(A[i][n+1],mul(mul(inv,p[i]),mul(a[v],h[mm])));
			Inc(A[i][v],mul(inv,dec(1,p[i])));
		}
	}
	for(int re i=1;i<=n;++i){
		int p=i;for(;!A[p][i];++p);
		if(p!=i)for(int re j=i;j<=n+1;++j)std::swap(A[i][j],A[p][j]);
		p=power(A[i][i],mod-2);
		for(int re j=i;j<=n+1;++j)Mul(A[i][j],p);
		for(int re j=i+1;j<=n;++j)if(int t=A[j][i])
		for(int re k=i;k<=n+1;++k)Dec(A[j][k],mul(t,A[i][k]));
	}
	for(int re i=n;i;--i){
		b[i]=dec(0,A[i][n+1]);
		for(int re j=i+1;j<=n;++j)Dec(b[i],mul(b[j],A[i][j]));
	}
}

signed main(){
#ifdef zxyoi
	freopen("island.in","r",stdin);
#endif
	init();scanf("%d",&n),m=c[n];
	assert(n<=50);
	for(int re i=1;i<=n;++i)scanf("%d",p+i);
	for(int re i=1;i<=n;++i){
		int k,u;scanf("%d",&k);
		while(k--){
			scanf("%d",&u);
			G[i].push_back(u);
		}
	}
	for(int re i=1;i<=n;++i)f[i][c[i]]=1;
	for(int re i=2;i<=n;++i)
	for(int re j=0;j<=c[i];++j)
	for(int re l=1;c[l]<=j;++l)
	Dec(f[i][j],mul(C[i-1][l],f[i-l][j-c[l]]));
	
	for(int re i=0;i<=m;++i){
		for(int re j=i;j<=m;++j)
		Inc(g[i],mul(f[n][j],C[j][i]));
		g[i]=dec(C[m][i],g[i]);
	}
	for(int re i=0;i<m;++i)h[i+1]=mul(mul(g[i+1],i+1),power(mul(m-i,g[i]),mod-2));
	memset(f,0,sizeof f);
	for(int re i=m;~i;--i)if(g[i])gauss(i+1,f[i+1],f[i]);
	cout<<dec(f[0][1],1)<<"\n";
	return 0;
}

代码（复杂度$O(n^4)$）：

#include<bits/stdc++.h>
#define ll long long
#define re register
#define cs const

using std::cerr;
using std::cout;

cs int mod=1e9+7;
inline int add(int a,int b){a+=b-mod;return a+(a>>31&mod);}
inline int dec(int a,int b){a-=b;return a+(a>>31&mod);}
inline int mul(int a,int b){ll r=(ll)a*b;return r>=mod?r%mod:r;}
inline int power(int a,int b,int res=1){
	for(;b;b>>=1,a=mul(a,a))(b&1)&&(res=mul(res,a));
	return res;
}
inline void Inc(int &a,int b){a+=b-mod;a+=a>>31&mod;}
inline void Dec(int &a,int b){a-=b;a+=a>>31&mod;}
inline void Mul(int &a,int b){a=mul(a,b);}

cs int N=55,M=N*N>>1;

int fac[M],ifac[M];
int c[N];
inline void init(){
	fac[0]=fac[1]=ifac[0]=ifac[1]=1;
	for(int re i=2;i<M;++i)fac[i]=mul(fac[i-1],i);
	ifac[M-1]=power(fac[M-1],mod-2);
	for(int re i=M-2;i;--i)ifac[i]=mul(ifac[i+1],i+1);
	for(int re i=1;i<N;++i)c[i]=i*(i-1)>>1;
}

inline int C(int n,int m){return n>=0&&m>=0&&n>=m?mul(fac[n],mul(ifac[m],ifac[n-m])):0;}
inline int invC(int n,int m){return n>=0&&m>=0&&n>=m?mul(ifac[n],mul(fac[m],fac[n-m])):0;}

int n,m;
std::vector<int> G[N];
int p[N],f[N][M],g[M],h[M];
int a[N][N<<1],e[N][M],coef[N][N];

void gauss(){
	for(int re i=1;i<=n;++i){
		int p=i;for(;!a[p][i];++p);
		if(p!=i)for(int re k=i;k<=2*n+1;++k)std::swap(a[i][k],a[p][k]);
		p=power(a[i][i],mod-2);
		for(int re k=i;k<=2*n+1;++k)Mul(a[i][k],p);
		for(int re j=i+1;j<=n;++j)if(p=a[j][i])
		for(int re k=i;k<=2*n+1;++k)Dec(a[j][k],mul(p,a[i][k]));
	}
	for(int re i=n;i;--i){
		for(int re j=n;j>i;--j)
		for(int re k=n+1;k<=2*n+1;++k)Dec(a[i][k],mul(a[i][j],coef[j][k-n]));
		for(int re k=1;k<=n+1;++k)coef[i][k]=a[i][k+n];
	}
}

signed main(){
#ifdef zxyoi
	freopen("island.in","r",stdin);
#endif
	init();scanf("%d",&n),m=c[n];
	assert(n<=50);
	for(int re i=1;i<=n;++i)scanf("%d",p+i);
	for(int re i=1;i<=n;++i){
		int k,u;scanf("%d",&k);
		while(k--){
			scanf("%d",&u);
			G[i].push_back(u);
		}
	}
	for(int re i=1;i<=n;++i)f[i][c[i]]=1;
	for(int re i=2;i<=n;++i)
	for(int re j=0;j<=c[i];++j)
	for(int re l=1;c[l]<=j;++l)
	Dec(f[i][j],mul(C(i-1,l),f[i-l][j-c[l]]));
	
	for(int re i=0;i<=m;++i){
		for(int re j=i;j<=m;++j)
		Inc(g[i],mul(f[n][j],C(j,i)));
	}
	for(int re i=1;i<=m;++i)h[i]=mul(g[i],invC(m,i));
	for(int re i=m;i;--i)Dec(h[i],h[i-1]);
	for(int re u=1;u<=n;++u){
		Inc(a[u][u],1);
		Inc(a[u][2*n+1],1);
		int inv=power(G[u].size(),mod-2);
		for(int re v:G[u])
		Inc(a[u][v+n],mul(inv,p[u])),
		Dec(a[u][v],mul(inv,1-p[u]+mod));
	}
	gauss();
	for(int re i=1;i<=n;++i)e[i][0]=0;
	for(int re i=1;i<=m;++i)
	for(int re j=1;j<=n;++j){
		for(int re k=1;k<=n;++k)Inc(e[j][i],mul(coef[j][k],e[k][i-1]));
		Inc(e[j][i],coef[j][n+1]);
	}
	int ans=0;
	for(int re i=1;i<=m;++i)Inc(ans,mul(e[1][i],h[i]));
	cout<<dec(ans,1)<<"\n";
	return 0;
}