分析
为了方便，设 $n<m<l$，$V=nml$。
考虑容斥，求出至少有 $i$ 个极大数的方案数 $c_i$。
那么 $c_i$ 等于选出 $i$ 个极大数的方案数、$i$ 个极大数所在的平面的并填数的方案数、其它位置填数的方案数之积。
首先考虑求选出 $i$ 个数的方案数 $f_i$，不难得到
$$
f_i=\prod_{j=0}^{i-1}(n-j)(m-j)(l-j)
$$
在求接下来两个方案数之前，先考虑算 $i$ 个极大数所在平面的并的方块数 $g_i$，显然有
$$
g_i=V-(n-i)(m-i)(l-i)
$$
现在考虑求 $i$ 个极大数所在平面的并的填数的方案数 $h_i$。显然每个数需要满足小于控制它的极大数中最小的那个，我们称这个极大数实际控制这个数。将所有极大数从大到小考虑，则最大的那个极大数实际控制 $g_i-g_{i-1}-1$ 个数（不包括自己），且有 $g_i-1$ 个数可以填，则方案数为 $A_{g_i-1}^{g_i-g_{i-1}-1}=\frac{(g_i-1)!}{g_{i-1}!}$，且填完后变为子问题 $h_{i-1}$。所以
$$
h_i=\prod_{j=1}^{i}\frac{(g_j-1)!}{g_{j-1}!}
$$
其他位置填数的方案数就很好算了，为 $(V-g_i)!$。
注意上面漏算了需要乘上选出 $g_i$ 个数的方案数。所以有
$$
\begin{aligned}
c_i&={V\choose g_i}f_ih_i(V-g_i)!\\
&=\frac{V!}{g_i!(V-g_i)!}f_i\prod_{j=1}^i\frac{(g_j-1)!}{g_{j-1}!}(V-g_i)!\\
&=V!f_i\prod_{j=1}^i(g_j-1)!\prod_{j=0}^{i}\frac{1}{g_j!}\\
&=V!f_i\prod_{j=1}^i\frac{1}{g_j}
\end{aligned}
$$
我们实际上要算的是至少有 $i$ 个极大的数的概率 $d_i$，直接除以 $V!$ 即可
$$
p_i=f_i\prod_{j=1}^i\frac{1}{g_j}
$$
线性求逆元即可 $\mathcal{O}(n)$ 求出所有 $p_i$。
最后再考虑求答案。因为有
$$
p_i=\sum_{j=i}^n{j\choose i}ans_j
$$
二项式反演得到
$$
ans_k=\sum_{i=k}^n(-1)^{i-k}{i\choose k}p_i
$$
总时间复杂度 $\mathcal{O}(n)$。
代码
// ===================================
// author: M_sea
// website: https://m-sea-blog.com/
// ===================================
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
int read() {
int X=0,w=1; char c=getchar();
while (c<'0'||c>'9') { if (c=='-') w=-1; c=getchar(); }
while (c>='0'&&c<='9') X=X*10+c-'0',c=getchar();
return X*w;
}
const int N=5000000+10;
const int mod=998244353;
int qpow(int a,int b) { int c=1;
for (;b;b>>=1,a=1ll*a*a%mod) if (b&1) c=1ll*c*a%mod;
return c;
}
int fac[N],ifac[N];
int C(int n,int m) {
if (n<m) return 0;
return 1ll*fac[n]*ifac[m]%mod*ifac[n-m]%mod;
}
void init(int n) {
fac[0]=1;
for (int i=1;i<=n;++i) fac[i]=1ll*fac[i-1]*i%mod;
ifac[n]=qpow(fac[n],mod-2);
for (int i=n;i;--i) ifac[i-1]=1ll*ifac[i]*i%mod;
}
int f[N],g[N],invg[N];
int main() { init(5e6);
int T=read();
while (T--) {
int n=read(),m=read(),l=read(),k=read(),V=1ll*n*m%mod*l%mod;
if (n>m) swap(n,m); if (m>l) swap(m,l); if (n>m) swap(n,m);
if (k>n) { puts("0"); continue; }
f[0]=1;
for (int i=1;i<=n;++i)
f[i]=1ll*f[i-1]*(n-i+1)%mod*(m-i+1)%mod*(l-i+1)%mod;
for (int i=1;i<=n;++i)
g[i]=(V-1ll*(n-i)*(m-i)%mod*(l-i)%mod+mod)%mod;
int mulg=1;
for (int i=1;i<=n;++i) mulg=1ll*mulg*g[i]%mod;
invg[n]=qpow(mulg,mod-2);
for (int i=n;i;--i) invg[i-1]=1ll*invg[i]*g[i]%mod;
int ans=0;
for (int i=k;i<=n;++i) {
int w=1ll*f[i]*invg[i]%mod;
if ((i-k)&1) ans=(ans-1ll*C(i,k)*w%mod+mod)%mod;
else ans=(ans+1ll*C(i,k)*w)%mod;
}
printf("%d\n",ans);
}
return 0;
}