僅僅是對 \(O(n^4)\) 做法的一個記錄。

簡要題意

有 \(N\) 座島嶼，初始時沒有邊。每座島嶼都有一個概率值 \(p_i\) 和一個大小為 \(s_i\) 友好列表 \(A_i\) 。

小 \(c\) 站在 \(1\) 號島嶼，依次執行以下操作:

\((1)\) 設現在在島嶼 \(x\)，有 \(p_x\) 的概率產生一條圖中尚未存在的隨機無向邊，不會產生自環。

\((2)\) 如果此時所有島嶼仍未聯通，她會在當前點的友好列表中，等概率隨機選擇一個，走到那座島嶼上。並把不滿意度增加 \(1\)，然後重複 \((1)\)。否則就結束這個過程

求她的期望不滿意度，對一個 \(10^9+7\)

取模.

\(n\le 50\).

解題思路

思考後不難發現，狀態只與當前的位置，以及邊數相關，那麼可以設狀態：

\(e(x,i)\) 表示現在位於點 \(x\)，圖還差 \(i\) 條邊連通，期望還要走多少步。

則最終答案是：

\[Ans = \sum_{i=1}^me(1, i)\cdot P(\text{恰好 }i\text{ 條邊連通的概率}) \]

不難得到 \(e\) 之間的轉移式：

\[e(x,i)=1+\frac{1}{s_x}\sum_{x\rightarrow y}p_x\cdot e(y,i-1)+(1-p_x)\cdot e(y,i) \]

如果分層高斯消元是 \(O(n^5)\)

，考慮設 \(e(x,i)=c(x)+\sum b(x,y)\cdot e(y,i-1)\) 來遞推解決。

帶入轉移式可得：

\[\begin{aligned} e(x,i)=1+\frac{1}{s_x}\sum_{x\rightarrow y}p_x\cdot e(y,i-1)+(1-p_x)\cdot (\sum_{z=1}^nc(y)+b(y,z)\cdot e(z,i-1))\\ e(x,i)=1+\frac{p_x}{s_x}\sum_{x\rightarrow y}e(y,i-1)+\frac{(1-p_x)}{s_x}\cdot \sum_{x\rightarrow y}(c(y)+\sum_{z=1}^nb(y,z)\cdot e(z,i-1))\\ e(x,i)=1+\frac{(1-p_x)}{s_x}\cdot \sum_{x\rightarrow y}^nc(y)+\frac{p_x}{s_x}\sum_{x\rightarrow y}e(y,i-1)+\frac{(1-p_x)}{s_x}\cdot \sum_{x\rightarrow y}\sum_{z=1}^nb(y,z)\cdot e(z,i-1)\\ \end{aligned} \]

於是可以得到兩個方程：

\[c(x)=1+\frac{(1-p_x)}{s_x}\cdot \sum_{x\rightarrow y}^nc(y) \] \[b(x,z)=\frac{p_x}{s_x}[x\rightarrow z]+\frac{1-p_x}{s_x}\sum_{x\rightarrow y}b(y,z) \]

這樣就可以在 \(O(n^4)\) 的高斯消元內得到所有的 \(b,c\)，然後遞推得到 \(e(1,i)\)。

---

接下來考慮如何計算恰好 \(i\) 條邊連通的概率，其實要算的就是 \(n\) 個點 \(i\) 條邊連通圖的方案數。

這類問題的經典處理方式是用斯特林反演來容斥。

於是考慮劃分連通塊，不同塊之間的邊不能選擇，同一塊內的邊可選可不選。

令至少 \(n\) 個連通塊的方案數為 \(f(n)\)，恰好 \(n\) 個連通塊的方案數為 \(g(n)\)，兩者關係即為：

\[\begin{aligned} f(n)&=\sum_{i=n}\begin{Bmatrix}n\\i\end{Bmatrix}g(i)\\ g(n)&=\sum_{i=n}(-1)^{i-n}\begin{bmatrix}i\\n\end{bmatrix}f(i) \end{aligned} \]

而我們希望通過計算 \(f\) 得到 \(g(1)\)，故 \(f(i)\) 的容斥係數為 \((-1)^{i-1}(i-1)!\)。

所以 \(i\) 條邊， \(n\) 個點的連通圖數量實際上為：

\[\sum_{G}(-1)^{c(G)-1}(c(G)-1)!{e(G) \choose i} \]

其中 \(G\) 需要取遍所有的 \(n\) 個點的，劃分成了若干連通塊的，每個連通塊內連成了完全圖的圖，\(c(G)\) 表示 \(G\) 的連通塊數量，\(e(G)\) 表示 \(G\) 的邊數。

設 \(dp(i,j)\) 表示 \(i\) 個點，劃分成了若干連通塊，每個連通塊內部連成完全圖後邊數和為 \(j\)，的所有圖帶上容斥係數的 \((-1)^{c(G)-1}(c(G)-1)!\) 之和，於是考慮 \(dp(i,j)\) 的轉移。

首先不難滿足 \((-1)^{c(G)-1}\) 這個係數條件，轉移的時候直接乘上 \(-1\) 的係數即可。

那麼 \((c(G)-1)!\) 怎麼辦呢，注意到這個係數的含義在於，我們希望把同一種劃分連通塊的方式計算 \((c(G)-1)!\) 次，那麼不妨賦初值時先確定 \(1\) 號點所在連通塊大小，轉移的時候列舉下一個連通塊的大小為 \(k\)，以 \(\dbinom{i+k-1}{k}\) 的係數重標號，這樣除了 \(1\) 號點所在連通塊位置一定排在第一個，其他的連通塊會因為在轉移時對於其新增順序有區分，有 \((c(G)-1)!\) 種排列順序，滿足了我們的需要。

於是有如下轉移式：

\[dp(i+k,j+\frac{k\cdot(k-1)}{2})\gets^{-} \sum_{k=1}^{n-i} dp(i,j)\cdot \binom{i+k-1}{k} \]

這樣可以 \(O(n^4)\) 計算出 \(dp(n,i)\) ，然後 \(O(n^4)\) 直接算出 \(g(i)\)，表示 \(n\) 個點 \(i\) 條邊的連通圖個數，那麼也不難算恰好 \(i\) 條邊連通的概率了。

namespace Gauss{
	int n, a[N][N], b[N];
	inline int get(int x){ return (LL) b[x] * qpow(a[x][x], mod - 2) % mod; }
	void init(int m){
		n = m;
		lfor(i, 1, n) b[i] = 0;
		lfor(i, 1, n) lfor(j, 1, n) a[i][j] = 0;
	}
	void solve(){
		lfor(i, 1, n){
			int inv = qpow(a[i][i], mod - 2);
			lfor(j, 1, n) if(i != j){
				int K = (LL) inv * a[j][i] % mod;
				lfor(k, 1, n) MOD(a[j][k] -= (LL) K * a[i][k] % mod);
				MOD(b[j] -= (LL) K * b[i] % mod);
			}
		}
	}
}
using Gauss :: get;
using Gauss :: init;
using Gauss :: solve;

void get_coef(){
	init(n);
	lfor(x, 1, n){
		Gauss :: b[x] = mod - 1;
		Gauss :: a[x][x] = mod - 1;
		int v = (LL) q[x] * qpow(s[x], mod - 2) % mod;
		lfor(j, 1, s[x]) Gauss :: a[x][a[x][j]] = v;
	}
	solve();
	lfor(x, 1, n) c[x] = get(x);

	lfor(z, 1, n){
		init(n);
		
		lfor(x, 1, n){
			int v = (LL) q[x] * qpow(s[x]) % mod;
			if(lnk[x][z]) Gauss :: b[x] = Mod(-(LL) p[x] * qpow(s[x]) % mod);
			Gauss :: a[x][x] = mod - 1;
			lfor(y, 1, n) if(lnk[x][y]) Gauss :: a[x][y] = v;
		}
		
		solve();
		lfor(x, 1, n) b[x][z] = get(x);
	}

	lfor(i, 1, m) lfor(x, 1, n){
		lfor(y, 1, n) MOD(e[x][i] += (LL) b[x][y] * e[y][i - 1] % mod - mod);
		MOD(e[x][i] += c[x] - mod);
	}
}

void get_prob(){
	lfor(i, 1, n) dp[i][C(i, 2)] = 1;
	lfor(i, 1, n) lfor(j, 0, m) if(dp[i][j]) lfor(k, 1, n - i) 
		MOD(dp[i + k][j + C(k, 2)] += -(LL) dp[i][j] * C(i + k - 1, k) % mod);
	lfor(i, 0, m) lfor(j, i, m)
		MOD(g[i] += (LL) dp[n][j] * C(j, i) % mod - mod);
	lfor(i, 1, m) g[i] = (LL) g[i] * qpow(C(m, i), mod - 2) % mod;
	rfor(i, m, 1) MOD(g[i] -= g[i - 1]);
}

signed main(){
	read(n), m = n * (n - 1) / 2;
	lfor(i, 1, n) read(p[i]), q[i] = Mod(1 - p[i]);
	lfor(i, 1, n){
		read(s[i]), a[i] = new int[s[i] + 1];
		lfor(j, 1, s[i]) read(a[i][j]), lnk[i][a[i][j]] = 1;
	}
	
	prep();
	get_coef();
	get_prob();
	
	lfor(i, 1, m) MOD(Ans += (LL) e[1][i] * g[i] % mod - mod);
	cout << Mod(Ans - 1) << endl;
	return 0;
}