來嘗試一下從容斥的最底層來解決這道題目，已經看過巨佬推過了，只是自己再推一次加深印象。

首先我們定義三元組 \(p(x,y,z)\) 表示一個點，再定義集合 \(S={p_1,p_2,...,p_k}\) 表示極大點的集合。

然後我們輕易得到對於這若干個極大點，不能出現兩兩存在相同的座標，這是一個限制。

易得，\(|S|_{\max}=\min(n,m,l)\) 。

然後我們定義 \(F_S\) 表示欽定集合 \(S\) 中的點是極大點，其他點隨便選的概率，我們考慮如何計算。

顯然的兩兩極大點之間的概率並不是獨立的，但是由於極大點是兩兩不等的（任意兩點都是不等的），所以我們可以欽定一個大小排列順序，我們按照從小到大的順序加入，由於交集的部分已經滿足了比小的小，所以我們不需要再欽定其大小關係了。

\[F_S=|S|!\prod_{i=1}^{|S|}\frac1{nml-(n-i)(m-i)(l-i)} \]

然後根據這個式子，我們可以輕易得到概率與具體的位置無關，只與 \(|S|\) 有關。

定於 \(f_i\) 為大小為 \(i\) 的 \(S\) 的 \(F_S\) 之和。

\[f_i=\sum_{|S|=i}F_S\\ =\sum_{|S|=i} i!\prod_{j=1}^{i}\frac1{nml-(n-j)(m-j)(l-j)}\\ =n^{\underline i}m^{\underline i}l^{\underline i}\prod_{j=1}^{i}\frac1{nml-(n-j)(m-j)(l-j)}\\ \]

我們設 \(g_i\)

表示恰好有 \(i\) 個極大點的概率，答案即是 \(g_k\) 。

\[f_i=\sum_{j=i}^{\min(n,m,l)}\binom{j}{i}g_j\\ g_i=\sum_{j=i}^{\min(n,m,l)}(-1)^{j-i}\binom{j}{i}f_j\\ \]

我們可以先將 \(\prod_{j=1}^{i}\frac1{nml-(n-j)(m-j)(l-j)}\\\) 預處理出來。然後再將 \(f_i\) 求出來，然後再將 \(g_k\) 求出來即可。

程式碼

#include<bits/stdc++.h>
using namespace std;
#define int long long
const int N=5e6+5;
const int MOD=998244353;
int n,m,l,k;
int fact[N],ifact[N],ksm_1[N];
int ksm(int x,int k){
	int res=1;
	for(;k;k>>=1,x=x*x%MOD)
	if(k&1) res=res*x%MOD;
	return res;
}
int dec_ksm(int x,int k){
	return fact[x]*ifact[x-k]%MOD;
}
int cal(int n,int m){
	if(n<0||m<0||n-m<0) return 0;
	return fact[n]*ifact[m]%MOD*ifact[n-m]%MOD;
}
int f[N],g[N],h[N];
void solve(){
	cin>>n>>m>>l>>k;
	int tot=min(n,min(m,l));
	g[0]=1;
	for(int i=1;i<=tot;++i) g[i]=g[i-1]*(n*m%MOD*l%MOD-(n-i)*(m-i)%MOD*(l-i)%MOD+MOD)%MOD;
	h[tot]=ksm(g[tot],MOD-2);
	for(int i=tot;i>=1;--i) h[i-1]=h[i]*(n*m%MOD*l%MOD-(n-i)*(m-i)%MOD*(l-i)%MOD+MOD)%MOD;
	f[0]=1;
	for(int i=1;i<=tot;++i) f[i]=h[i]*g[i-1]%MOD*f[i-1]%MOD;
	for(int i=1;i<=tot;++i) (f[i]*=dec_ksm(n,i)*dec_ksm(m,i)%MOD*dec_ksm(l,i)%MOD)%=MOD;
	int res=0;
	for(int i=k;i<=tot;++i) (res+=ksm_1[i-k]*cal(i,k)%MOD*f[i]%MOD)%=MOD;
	return printf("%lld\n",res),void();
}
signed main(){
	ksm_1[0]=1;
	for(int i=1;i<N;++i) ksm_1[i]=ksm_1[i-1]*(MOD-1)%MOD;
	fact[0]=1;
	for(int i=1;i<N;++i) fact[i]=fact[i-1]*i%MOD;
	ifact[N-1]=ksm(fact[N-1],MOD-2);
	for(int i=N-1;i>=1;--i) ifact[i-1]=ifact[i]*i%MOD;
	int T;cin>>T;
	while(T--) solve();
	return 0;
}