来尝试一下从容斥的最底层来解决这道题目，已经看过巨佬推过了，只是自己再推一次加深印象。
首先我们定义三元组 (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)}inom{j}{i}g_j\
g_i=sum_{j=i}^{min(n,m,l)}(-1)^{j-i}inom{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
",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;
}